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dynamic response measurements showed surprisingly agitated behavior, especially around 
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used to develop dynamic models to describe harmonic-drive operation.. From these 
models, I realized that non-linear frictional effects cannot be ignored in any accurate 
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kinematic error and compliant behavior are essential. By introducing the effects of friction, 
compliance, and kinematic error into a simplified representation of the gear-tooth meshing 
geometry in the actual transmission, model performance, especially around system 
resonance, improved substantially. Unfortunately, since accurate measurements of 
harmonic-drive properties were elusive and highly dependent on operating conditions, it 
seems unlikely that any harmonic-drive representation will be able to deliver reliable 
performance unless detailed experimental observations of the specific harmonic-drive 
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Chapter 1: Introduction 


A good mechanical designer holds a clear awareness of the capabilities and 
limitations of his design. For most designers, this awareness rests heavily on a firm 
understanding of the individual components that comprise the entire assembly. For 
mechanical systems employing harmonic drives, a designer can gain confidence in the 
performance of the overall mechanism if the operation of the harmonic drive is clearly 
understood. Additionally, for applications, such as robot control, which require 
mathematical models to describe system performance, a thorough understanding of the 
behavior of harmonic-drive transmissions is essential. The purpose of my investigation is 
to identify important attributes that influence the operation of harmonic-drive transmissions 
and to relate these attributes to overall dynamic behavior. On the basis of insight gained 
from this empirical analysis, I hope to develop mathematical models that can be used to 
predict reliably the dynamic performance of harmonic-drive transmissions. 

I began my investigation of harmonic-drives with an experimental analysis of 
kinematic, static, and dynamic behavior. In particular, I recorded and analyzed kinematic 
position error, compliance, and frictional losses in the transmission. From these results I 
discovered that friction and flexibility in the harmonic drives can be highly non-linear and 
very sensitive to changes in environmental and operating conditions. Because of these 
poorly behaved transmission properties, the dynamic response of the transmission was 
often equally erratic and unpredictable. Specifically, operation in regions that excited 
system resonance bred significant energy losses and surprisingly turbulent behavior. 
Despite the complexity of these experimental observations, I undertook to develop simple 
harmonic-drive representations to describe the discovered behavior. In total, five different 
models were derived and evaluated to illustrate the range of accuracy that could be delivered 
by harmonic-drive representations of varied complexity. From these models, the 
importance of friction, compliance, and kinematic-error on model performance was 
emphasized as was the importance of gear-tooth meshing geometry in the transmission. 
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Chapter 1: Introduction 


Results from these models were encouraging and informative, but more detailed 
representations are required to capture fully the subtleties of harmonic-drive behavior. 

In this document, I will describe completely the results of my experimental and 
theoretical analysis of harmonic drives. I will begin this presentation with a summary of 
existing research that is related to my investigation. Second, I will display and discuss the 
results of my experimental inquiry into important transmission properties and dynamic 
response. Third, I will develop several theoretical models to describe the observed 
behavior, and last, I will offer conclusions and recommendations to guide future work. 


1-1 The Harmonic Drive 

Since its conception by Walt Musser in 1955, the harmonic drive has found 
widespread use and acceptance among mechanical designers. This nascent mechanical 
transmission, occasionally labeled “strain-wave gearing”, employs a continuous deflection 
wave along a non-rigid gear to allow gradual engagement of gear teeth. Because of this 
unconventional gear-tooth meshing action, harmonic drives can deliver very high reduction 
ratios in a very small package. As a peripheral consequence, the radical mechanical 
operation of this gear train defies conventional understanding of gear behavior and creates a 
new arena for exploration and understanding. Following this impetus for investigation, 
several independent researchers have made valuable contributions to the state of harmonic- 
drive technology in the past thirty years. In addition to surveying the operational principles 
and utility of harmonic transmissions, this section will attempt to present a brief yet broad 
summary of existing research insights related to my investigation. 


1.1.1 The Principles of Operation 

Every harmonic drive is composed of the three distinct parts illustrated in figure 
1.1: the wave-generator, the flexspline, and the circular spline. The wave-generator is a 
ball-bearing assembly with a rigid, elliptical inner-race and a flexible, thin-walled outer 
race. The flexspline is a thin-walled, flexible cup adorned with small, external gear teeth 
around its rim. When assembled, the wave-generator is nested inside the flexspline, 
causing the flexible gear-tooth circumference on the flexspline to adopt the elliptical profile 
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wave-generator ■ flexspline ■ circular spline 


Figure 1.1: The three components of a harmonic drive 



Figure 1.2: An assembled harmonic-drive transmission 
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of the wave-generator. The circular spline is a rigid ring with internal teeth machined along 
a slightly larger pitch diameter than those of the flexspline. When the flexspline-and-wave- 
generator assembly is inserted into the circular spline, the external teeth on the flexspline 
mesh with the internal teeth on the circular spline along the major axis of the wave- 

generator ellipse. A correctly assembled harmonic drive with a cut-away section is shown 
in figure 1.2. 

If properly assembled, all three components of the transmission can rotate at 
different but coupled velocities on the same axis. The rotational mechanism, as illustrated 
in figure 1.3, begins when rotation of the wave-generator carries the zone of gear-tooth 
engagement with its major axis. When this zone is propagated 180° around the 
circumference of the circular spline, the flexspline, which contains two fewer teeth than the 
circular spline, will lag by one tooth relative to the circular spline. Through this gradual 
and continuous engagement of slightly offset gear teeth, every rotation of the wave- 
generator moves the flexspline two teeth back on the circular spline. Through this 
unconventional mechanism, gear ratios up to 320:1 can be achieved in a single 
transmission. 


Since the harmonic drive has three rotational pons, it boasts a versatility unavailable 
to standard gear transmissions. Specifically, in its most popular configuration, the circular 
spline is fixed to ground and a low-torque, high-speed motor driving the wave-generator 
can produce high-torque, low-speed rotation on the flexspline. Similarly, with the 
flexspline mounted to ground, the torque on the wave-generator can be magnified and 
transmitted through the circular spline. If either the circular spline or flexspline is used as 
an input, the wave-generator can be driven at high velocities and low torques. In general, 
by using different combinations of rotations on the three harmonic-drive components, 
numerous differential-gearing functions and reduction ratios can be achieved. 


1.1.2 Performance Features and Applicability 

Currently, commercial harmonic drives are manufactured by two companies: The 
Harmonic Drive Machinery Group in Wakefield, MA, USA, and Harmonic Drive Systems, 
Incorporated, in Tokyo, Japan. Available harmonic drive products include the standard 
cup-type gearing, pancake transmissions, housed units, and motor-driven motion-control 
systems. Typical gear ratios of commercial drives range from 50:1 to 320:1, and 
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Figure 1.3: The harmonic-drive operating mechanism 
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efficiencies can approach 90 percent. The smallest transmissions can provide a respectable 
maximum torque output of about 3.5 N-m (30.0 in-lbs), while the heavy-duty units boast 
up to 10,000 N-m (89,000 in-lbs) of torque capacity. 

Because of its unique operating principles, the harmonic drive displays performance 
features both superior and inferior to conventional gear transmissions. These advantages 
and disadvantages are itemized in Table 1.1 and Table 1.2. 


The unique performance features of the harmonic drive have captured the attention 
of designers in many fields. Specifically, due to their accuracy and simple construction. 


Table 1.1: Performance advantages of harmonic. drive transmissions 

High torque capacity 

Since torque is transmitted through multiple-tooth contact, harmonic 
drives can withstand high loads at small pitch-diameters. 

Concentric geometry 

Since all three harmonic-drive components are concentric and coaxial, 
designers can drastically reduce power-train size and complexity. 

Lightweight and 
compact design 

Requiring only three basic elements, the harmonic drive can deliver 
extremely high gear ratios in a small package. 

Zero backlash 

Natural gear preloading and predominantly radial tooth-engagement 
eliminate virtually all transmission backlash. 

High Efficiency 

If properly lubricated, typical efficiencies of harmonic-drive 
transmissions can reach 80 to 90 percent. 

Backdrivability 

Due to their high efficiency, wave-generator rotation can be driven 
through the flexspline or circular spline. 


Table 1.2: Performance disadvantages of harmonic-drive transmissions 


High flexibility Due to the high loads seen by the wave-generator and gear teeth, 

moderate operating torques can produce substantial transmission 
torsion. 

Kinematic error Due to manufacturing inaccuracies, harmonic drives exhibit a small but 

ubiquitous position error across the transmission. 

Resonance vibration Since torque fluctuations produced by kinematic error can interact with 

the low stiffness of the transmission to excite resonances, high 
vibration amplitudes may be generated in some operating ranges. 

Non linearity Both the flexibility and frictional losses in the drive exhibit hiqhly 

nonlinear behavior. 

Poorly understood Compared to conventional gear transmissions, relatively little 

exploration has been done of the unusual operating mechanisms in 
the harmonic drive. 
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harmonic drives have been used successfully in industrial robots, assembly equipment, 
and measuring instruments. Heavy-duty applications such as machine tools and printing 
presses have also utilized harmonic drives for their high-torque capabilities. Additionally, 
space and aircraft systems often employ harmonic drives for their lightweight and compact 
geometry. 


1.1.3 The State of Harmonic-Drive Technology 

Throughout its short existence, the harmonic drive has enjoyed continuously 
increasing international attention from designers as well as researchers. In the Soviet 
Union, substantial research, predominantly performed on a harmonic drive featuring a 
wave-generator with two rollers instead of an elliptical bearing, has provided valuable 
insight into the structural and dynamic properties of the drive as well as its suitability for 
heavy-duty applications. More recently, Japanese researchers have explored such areas as 
transmission stiffness, positioning accuracy, and tooth-meshing mechanisms in an effort to 
develop new gear-tooth geometries and improve harmonic-drive performance. 
Investigating harmonic drive performance in robotics and aerospace applications, 
researchers in the United States have also contributed constructive insights about the limits 
of harmonic-drive applicability. 

The reference list at the end of this document contains citations for numerous papers 
encompassing the scope harmonic-drive research worldwide. Primarily, this research has 
been focused in six major areas: (1) vibration and dynamic behavior, (2) kinematic error, 
(3) torsional stiffness, (4) structural analysis, (5) robotics applications, and (6) gear-tooth 
geometry. Since my investigation is particularly focused on the dynamic and static 
performance of harmonic drives, past research on dynamic behavior, positional error, and 
torsional stiffness received most of my attention. In an effort to overview the current state 
of harmonic drive technology, a description of the fundamental mechanisms of harmonic- 
drive operation as well as relevant insights presented in existing literature will be 
summarized in this section. 


1.1.3.1 Fundamental Operating Principles 

In commercial catalogs, every harmonic drive is assigned a transmission ratio, N, 
which describes its position, velocity, and torque behavior. Specifically, given a known 
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rotation on two of the three harmonic-drive ports as well as a value for N, the ideal rotation 
of the third harmonic-drive port can be predicted by the equation 

0 wg = (N + 1) 0 CS - N 0 fs , (i.i) 

where 0 wg is the rotation of the wave-generator, 0 CS is the rotation of the circular spline, 
and 0f s is the rotation of the flexspline. All three rotations in this equation are defined in 

the same frame of reference. Similarly, given that N is constant under ideal assumptions, 
the derivative of this relationship yields a similar velocity constraint: 

co wg = (N + 1) co cs - N C0f s , (12) 

where co wgf co CSt and C 0 f s represent the angular velocities of the three harmonic-drive 
components relative to the same velocity reference. From these equations, it can be seen 
that, if the circular-spline velocity is zero, the wave-generator will rotate N-times faster than 
the flexspline in the opposite direction. In a different configuration, a grounded flexspline 
dictates that the wave-generator will spin (N+l) times faster than the circular spline in the 
same direction. By applying the law of power conservation to the three ports of the 
harmonic drive, the ideal torque behavior of the transmission can also be defined: 

Twg = (N J Mj Tcs = 'N Tfs> d- 3 ) 

where T wg , T cs , and T fs are the torques, defined with identical sign conventions, seen by 
the three harmonic-drive components. From this identity, it can be seen that the torque on 
one port of the harmonic drive dictates the ideal torque seen by the other two ports. Notice 
that the torque on the flexspline is nearly equivalent to the torque on the circular spline, 
which is approximately N-times greater than the torque on the wave-generator. By 
applying the torque, velocity, and position equations presented above, the dynamic 

behavior of a harmonic drive transmission under ideal assumptions can be completely 
defined. 


1.1.3.2 Kinematic Error 

In applications requiring high positional accuracy, the inherent kinematic errors 
manifested by harmonic drives expose the shortcomings of an ideal transmission model. 
This position error, 0 err , is typically measured by subtracting the rotation at the output of 
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the harmonic drive from the input rotation scaled by the ideal gear ratio for the given 
transmission configuration: 

ft - ®in _ Q 
u err — • t'out • 

gear ratio (\ 4 ) 

Based on experimental observation, manufacturers’ catalogs report that the magnitude of 
this position error typically varies periodically at a frequency of twice the wave-generator 
rotational velocity. Additionally, the amplitude of this periodic error often modulates twice 
per flexspline or circular-spline revolution as illustrated in figure 1.4. One manufacturer, 

[3] , reports that the peak-to-peak amplitude of the kinematic error on all harmonic drive 
models rarely exceeds 0.033 degrees of output rotation (2 arc minutes), while the other, 

[4] , claims that peak-to-peak error is inversely proportional to pitch-diameter and ranges 
from 0.095 degrees (5.7 arc minutes) to 0.013 degrees (0.8 arc minutes). Both 
manufacturers claim that position-error can be reduced further on a custom basis. [ 3 , 4 ] 

The root of harmonic-drive kinematic inaccuracy has been attributed to several 
different sources. As reported by Nye and Kraml in [24], harmonic-drive manufacturers 
claim that the position error is caused predominantly by manufacturing and assembly 
imperfections. Specifically, tooth-placement errors on both the circular spline and 
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flexspline, out-of-roundness in the three transmission components, as well as misalignment 
during assembly can account for the majority of kinematic imperfection. By developing a 
detailed model which translated expected manufacturing and assembly errors into overall 
transmission error, Emel’yanov, et al., [20], made accurate predictions of experimentally 
observed position error. Klypin, et al., [23], used Emal’yanov’s model to develop a 
computer simulation and further confirm the presence of numerous manufacturing and 
assembly errors in the measured transmission error. More recently, Yanabe, et al., [28], 
used precise experimental equipment to measure tooth-placement errors and concluded that 
the origins of position error lie mostly in the angular pitch errors of the flexspline and 
circular spline. In contrast to these conclusions, investigations by Hsia, [21, 22], and 
Ahmadian, [19], have focused blame for kinematic inaccuracies on inherent errors in the 
operating principles of harmonic drives irrespective of manufacturing and assembly errors. 
Both researchers developed kinematic models to predict position error based on 
transmission geometry and operating behavior. 

In addition to postulating probable sources for kinematic error, the experimental 
investigations of Yanabe, et al., [28], and Nye and Kraml, [24], have lent further insight 
into the repeatability and variation of typical gear-error waveforms. By measuring the 
position error of multiple harmonic drives over one complete output rotation, both 
researchers observed that the resulting error signature can display frequency components at 
two cycles per wave-generator revolution and several subsequent harmonics. Typically, 
error components occurring once per wave-generator revolution are not found since they 
are canceled by the symmetry of the wave-generator, which ensures two-point contact 
between of the flexspline and circular spline. As Nye and Kraml summarized, the primary 
error component was seen at two cycles per wave-generator revolution with the next nine 
harmonics contributing up to fifty-five percent of the total kinematic error. The authors 
regret that there is no method for predicting the frequency components of the position error 
signal for an arbitrary harmonic drive other than careful measurement. By exposing a 
range of harmonic drive models to a variety of environmental and operational conditions, 
Nye and Kraml catalogued variations in position-error due to (1) temperature, (2) wear, (3) 
pitch diameter, (4) misalignment, (5) improper assembly, and (6) overloading. Based on 
this new knowledge of gear-error behavior, Nye and Kraml then successfully developed an 

electronic motion-control scheme designed to compensate for known kinematic 
inaccuracies. 

Additional observations in [28] and [24] record the presence of the significant 
amplitude modulation or beating in position error waveforms of most harmonic drives 
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while a few transmissions showed little or no amplitude modulation. Both authors concur 
that the presence of amplitude modulation over several wave-generator cycles is caused by 
the summation of the circular spline and flexspline gear-error components; since these two 
error waveforms appear at slightly different frequencies, constructive and deconstructive 
interference of the two error signals is manifested by low-frequency beating. 


1.1.3.3 Torsional Stiffness 

Under ideal assumptions, a harmonic-drive transmission is treated as a perfectly 
rigid gear reduction. However, due to the relatively low torsional stiffness of harmonic 
drives, a more detailed understanding of transmission flexibility is often required for 
accurate modeling. Available research results have focused on understanding the 
characteristic shape and variation of the non-linear stiffness profile and pinpointing sources 
of transmission flexibility. 


Harmonic drive stiffness is commonly measured by locking the wave-generator and 
measuring the rotation of the output (circular spline or flexspline) as the result of an applied 
load. As described in a manufacturer’s catalog, [3], the typical shape of the stiffness curve 
for a harmonic drive loaded in both direction is shown in figure 1.5. This curve illustrates 
two characteristic properties of harmonic-drive flexibility: increasing stiffness with 



Figure 1.5: Typical stiffness profile displayed by harmonic drives 
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displacement and hysteresis loss. To capture the non-linear stiffness behavior, 
manufacturers, [3, 4], suggest using piecewise linear approximations whereas several 
independent researchers, [11, 18], prefer a cubic polynomial approximation. For certain 
modeling applications, other researchers, [42, 52] have had success using a linear stiffness 
approximation. At low applied torques, soft wind-up in the transmission due to this non¬ 
linear stiffness behavior can sometimes act like backlash. Fortunately, manufacturers claim 
that this deceptive non-linearity is kept below 0.05 degrees of output rotation (3 arc 
minutes). The hysteresis loss in a harmonic drive is a phenomenon more difficult to model 
than the stiffness, and consequently, it is ignored often. In support of this neglect, the 
manufacturer’s catalog, [3], reassures modelers that the total loss is typically less than 
0.033 degrees of output rotation (2 arc minutes). 

The primary sources of flexibility in harmonic-dnve transmissions has been a topic 
of much investigation and speculation. By performing tests on high-torque harmonic 
drives with two-cam wave-generators, Margulis, et al., [29], concluded that the largest 
component in the total torsional compliance was the deformation of the wave-generator. 
Specifically, it was observed that high radial forces on the wave-generator can cause 
substantial radial deflection allowing relative rotation across the transmission. By 
performing similar stiffness tests on conventional harmonic drives with elliptical-bearing 
wave-generators, Yanabe, et al., [32], also noticed that, due to the high radial forces seen 
by the wave-generator, resulting deflection in the wave-generator bearing can greatly 
influence transmission torsion. The zones of maximum loading on the wave-generator 
were predicted using a finite-element model of the flexspline. Research by Nye, [65], 
attributed the characteristic non-linear profile of harmonic-drive stiffness curves to 
increasing gear-tooth contact area under applied loads. Qualitative observations indicated 
that the axial length of the contact interface between meshing gear teeth can vary from ten 
percent, under no-load conditions, up to 100 percent, when fully loaded. The importance 
of gear-tooth contact on transmission stiffness was confirmed by Kiyosawa, et al., [57], 
when flexibility was dramatically reduced in harmonic drives with a modified gear-tooth 
profile designed to increase the tooth-engagement zone. In summary, existing research 
indicates that most transmission compliance can be traced to wave-generator deformation, 
but gear-tooth contact-area and meshing-zone, which can effect the load-distribution inside 
the harmonic drive, are influential in determining the shape of the stiffness profile. 

As a consequence of the relatively complex compliance mechanisms of the 
harmonic drive, significant variations can occur in transmission stiffness for a given drive. 
In particular, a manufacturer’s catalog, [3], cautions that stiffness profiles measured from 
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a single drive can vary as much as ±30 percent. Nye, [65], confirms the magnitude of this 
variability and observes its dependence on the integrity of internal fit-up. Specifically, he 
notes that stiffness typically increases with tighter meshing of gear teeth and increased 
internal loading. Experimental results from Nye and Hikada, et al., [11], indicate that, 
since wave-generator orientation often influences gear-tooth loading, transmission stiffness 
can vary significantly with wave-generator rotation. Further investigation by Nye 
produced theoretical loading models which successfully described the influence of gear- 
tooth preloading on stiffness and predicted preload- and stiffness-relaxation due to gear- 
tooth wear. 


1-1.3.4 Frictional Losses 

Another performance metric in which harmonic-drive behavior falls short of ideal is 
frictional dissipation. Manufacturer’s catalogs, [3, 5], report typical transmission power 
losses greater than ten percent. This transmission inefficiency is manifested by lower 
output torques than ideal assumptions would predict. Further research by Mariler and 
Richard, [52], among others, reports that the velocity-dependent damping in harmonic- 
drive transmissions increases less sharply at higher velocities than lower velocities and 
might even decrease at very high velocities. Based on these observations, the authors 
recommend using a cubic approximation to capture this non-linear damping behavior. 
Research has also been performed to understand the behavior of velocity-independent 
losses, such as static and Coulomb friction, in harmonic drives. Specifically, Nye, [65], 
reports that starting torque values for harmonic-drive transmission can fluctuate greatly as a 
function of the orientation of the wave-generator. The author attributes this variation to 
different levels of gear-tooth preloading causing fluctuating static friction torques over one 
revolution of the wave-generator. Static and kinematic models developed by Nye and 
others, [39], rely on geometric models of gear-tooth rubbing to estimate Coulomb friction 
forces on the gear teeth and wave-generator. 


1.1.3.5 Dynamic Behavior 

Due to the non-ideal behavior of harmonic-drive kinematic accuracy, stiffness, and 
frictional losses, ideal-transmission models are often insufficient to describe accurately the 
dynamic behavior of harmonic drives. In the past, many researchers, such as [42], have 
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opted to capture harmonic-drive dynamics using linear stiffness and damping relationships 
while others, [44, 51, 52], required non-linear representations to ensure model accuracy. 
In order to capture the torque fluctuations induced by kinematic error, a few researchers, 
such as [11, 18], supplement their models with a sinusoidal torque function derived from 
position-error fluctuations or measured directly. For many applications in which non-ideal 
transmission behavior becomes dominant, harmonic-drive manufacturer’s frequently 
recommend adjustment of the operating envelope to a region less plagued by transmission 
dynamics. Some research, such as [13], has been focused on developing simple models to 
predict operating regions in which harmonic drives are well-behaved. 

Depending on the design and operating conditions, several researchers note that 
strikingly non-ideal behavior can be exhibited in harmonic-drive systems. Specifically, 
experimental measurements by Hikada, et al., [11], and Volkov, et al, [18], reveal that 
harmonic drives typically exhibit vibrations at twice the frequency of wave-generator 
rotation and several subsequent harmonics. It was also noticed that the amplitude of 
vibration increased when one of these frequencies coincided with the natural frequency of 
the harmonic-drive system. Furthermore, plots of resonance vibration presented by 
Hikada, et al., illustrated the same twice-per-output-revolution amplitude modulation seen 
in position-error signatures. By developing a simple non-linear vibration model which 
captured the behavior of tooth-meshing errors, Hikada, et al., illustrated the direct influence 
of position error on dynamic vibration. Further comparison between model predictions and 
experimental results revealed the importance of representing the variation in stiffness with 
wave-generator angle to ensure model accuracy. Another non-linear vibration model 
developed by Volkov, et al., was able to predict observed system non-linearities such as 
vibrations at subharmonics of the wave-generator excitation frequency due to the non-linear 
harmonic-drive stiffness. 


1.2 The Scope of My Investigation 

The results of my research yielded several insights which may be useful to 
harmonic-drive users and designers. From my experimental analysis, I characterized the 
magnitude and typical behavior of several important transmission properties. In particular, 
kinematic error was confirmed to appear primarily at frequency-multiples of the wave- 
generator angular velocity, and static stiffness measurements were found to yield non-linear 
profiles that were corrupted by transmission friction. Additional measurements of static 
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and dynamic friction demonstrated that energy dissipation in harmonic drives varies non- 
linearly with velocity and periodically with rotation and can escalate during system 
resonance. By comparing these empirical observations to manufacturer’s estimates, I 
discovered that, due to the fickle nature of most harmonic-drive attributes, catalog 
predictions often fail to describe the observed behavior. Following this investigation of 
individual properties, overall dynamic response measurement showed surprisingly agitated 
transmission behavior. Specifically, due to the problematic influence of system resonance 
in typical harmonic-drive operating ranges, dynamic response data demonstrated violent 
torque and velocity fluctuations as well as sudden unexpected increases in operating 
velocity. Using the understanding of transmission properties gained from experimental 
analysis, this unusual dynamic behavior was rationalized. 

Using the experimental findings as a reference. Five different theoretical models 
were developed to describe harmonic-drive operation. First, an ideal transmission 
representation was constructed that showed very poor performance under dynamic 
conditions. When non-linear frictional losses were appended to this ideal representation, 
model integrity improved. In order capture the important influence of resonant vibration on 
system performance, non-linear compliance and kinematic error were introduced into the 
harmonic-drive representation. The resulting model demonstrated resonant frequencies and 
resulting vibrations close to experimental observations but remained ignorant to the 
enhanced frictional losses at resonance. By developing a simplified geometric 
representation of gear-tooth meshing that included Coulomb friction at the gear-tooth 
interface, the elusive resonance losses were replicated and model performance improved. 
Unfortunately, due inaccuracies that remained in this model and the difficulty of estimating 
reliable values for transmission properties, a more detailed representation of gear-tooth 
meshing behavior is required to improve model integrity. Nevertheless, the relatively 
simple mathematics of this representation and its close conceptual relationship to the actual 
transmission makes it a valuable tool for fostering an intuitive understanding of harmonic- 
drive operation. 


1.3 The Organization of this Document 

The presentation of research findings in this document begins with a complete 
discussion of my experimental investigation. First, I will describe the testing apparatus and 
instrumentation followed by a discussion of four important harmonic-drive properties: (1) 
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kinematic-error, (2) stiffness, (3) starting torque, and (4) dynamic friction. For each of 
these topics, the testing procedure will be described, experimental results will be presented, 
and conclusions will be made about the observed behavior. Using this empirical insight, 
dynamic-response measurements will then be characterized and rationalized in terms of 
important transmission properties. With experimental dynamic measurements to serve as a 
metric for model performance, my theoretical exploration of harmonic-drive operation will 
then be outlined. This presentation begins with a discussion of the dynamic models 
developed to simulate the testing apparatus followed by a description of five different 
harmonic-drive representations. For each of these transmission models, equations of 
motion will be derived, accuracy and reliability will be evaluated, and conclusions will be 
drawn about the potential usefulness and applicability of each representation. Based on 
results from these theoretical models as well as the findings of the experimental inquiry, 
conclusions and recommendations will be made about rewarding areas for future research. 
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A sound understanding of harmonic-drive behavior is vital for making well- 
informed design decisions about applications employing harmonic-drives. Additionally, 
the development of a reliable transmission model rests heavily upon accurate 
characterization of harmonic-drive performance parameters. Given this impetus for 
improved understanding, several tests have been performed to investigate the kinematic, 
static, and dynamic properties of this unique transmission. 

From my experimental investigation, I have collected data to characterize the 
kinematic error, stiffness, and friction behavior of three harmonic-dnve test specimens and 
have identified relationships between these transmission properties and overall harmonic- 
drive dynamic response. In general, I observed that, due to the inseparable interaction 
between friction, compliance, and transmission error, accurate characterization of these 
properties was difficult in many cases. Not surprisingly, because of this complication, 
catalog estimates of these values were often unreliable. Nevertheless, given the collected 
data, simple mathematical representations were constructed to extract useful information for 
dynamic modeling. Perhaps more importantly, from the qualitative understanding gained 
from the experimental analysis of harmonic-drive position error, stiffness, and friction, 
typicai harmonic-dnve dynamic response patterns were identified and rationalized in terms 
of the observed transmission properties. 

This section is dedicated to presenting the findings of my experimental investigation 
into these areas. To accomplish this task, the experimental apparatus used for all harmonic- 
dnve testing will be discussed first followed by a presentation of results from each of the 
experimental tests. During this presentation of the individual experiments, specific testing 
procedures will be outlined, results will be discussed and conclusions will be formed. The 
discussion of results in each of these sub-sections not only focuses on observed 
experimental behavior but also addresses experimental accuracy, catalog predictions, and 
potential modeling strategies. Finally, with the insight gained from all of the experimental 
tests, typical dynamic-response behavior observed on harmonic drives will be presented 
and characterized. 
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2.1 Description of Testing Apparatus 

All experimental testing for this project was performed on a three-joint, planar robot 
with harmonic-drive transmissions. This testing facility is shown in figure 2.1. To allow 
for straightforward data collection on the three harmonic drives, the apparatus was 
designed so that each joint could be removed from the robot and operated independently. 
For the purpose of this document, these three joints, each sporting a different harmonic 
drive, will be referenced by numbers 1, 2 and 3. This section is devoted to describing the 
design and instrumentation of these three joints followed by an overview of the data- 
acquisition methods used for each testing assembly. 


Hardware Design and Sensor Locations 


Due to the functional requirements of the robot, each of the three joint assemblies 
requires a different design. Serving as the robot base, joint 1 is the largest and most 
powerful of the three testing facilities. Joint 2 places second in bulk and brawn, relying on 
a slightly smaller motor and transmission than joint 1. With a much lower torque capacity 
than its two companions, joint 3 is significantly smaller in size than joint 1 or 2. Tables 2.1 
and 2.2 summarize important motor and transmission properties for each of the three joints, 

and in light of these differences, the specific design of these three testing stations can now 
be discussed. 


| Table 2.1: Specifications for the three harmonic drives tested 


4omt i 

Joint 2 

Joint a 

Model Ntofcer. * 

1M 

5C 

1C 

Catalog Trartsrofsaloo Ratio 

160:1 

160:1 

80:1 

Pftcft Diameter 

8.14 cm 

6.35 cm 

3.56 cm 


(3.20 in) 

--- 

(2.50 in) 

(1.40 in) 
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Table 2.2: Motor specifications for the three testing stations 




illilii 








Aerotech 

1035DC 

Aerotech 

1017DC 

Clifton SmCo 
AS-780 

1.84 N-m 
(16.25 in-lb) 

0.92 N-m 
(8.13 in-lb) 

0.35 N-m 
(3.13 in-lb) 

6000 rpm 

6500 rpm 

8000 rpm 

0.06 N-m/amp 
(0.53 in-lb/amp) 

0.03 N-m/amp 
(0.26 in-lb/amp) 

0.04 N-m/amp 
(0.36 in-lb/amp) 

6.30 volts/krpm 

3.04 volts/krpm 

4.25 volts/krpm 



Figure 2.2: Cross-section of joint 1 
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The design of the first harmonic-drive testing station, joint 1, is shown in figure 

2.2. From this cross-section, it can be seen that the flexspline of the harmonic drive is 
mounted to ground while the wave-generator is driven by a DC motor mounted to the 
circular spline. The output link of this joint is supported by a pair of angular-contact 
bearings and is driven by the circular spline of the harmonic drive. This joint also contains 
two encoders, one torque sensor, and a tachometer whose locations are illustrated in figure 

2.3. From this diagram, it can be seen that the input encoder and tachometer measure the 
rotation of the motor shaft with respect to the motor housing which is mounted to the 
circular spline. The output encoder is mounted to ground and measures the absolute 
rotation of the circular spline through a gear pair with an 8:1 reduction ratio. Lastly, since 
the only rigid connection between the drive train and ground is the torque sensor at the 
flexspline output, this sensor can be used to measure directly the torque across the 
transmission as seen by the flexspline. 
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As illustrated in Figure 2.4, the design of joint 2 is almost identical to joint 1. 
Specifically, on the input side of the joint, the wave-generator is driven by a motor 
mounted to the circular spline, while the rotation of the output link is driven by the circular 
spline. As seen on joint 1, the flexspline is attached to ground and angular-contact bearings 
are used to support the rotation of the output link. The instrumentation similarities between 
joint 1 and 2 are also apparent from the schematic in figure 2.5. As explained for joint 1, 
the input encoder and tachometer monitor the rotation of the motor shaft with respect to the 
circular spline while the output torque sensor measures all torque seen at the flexspline. 
Due to design restrictions, the output rotation on joint 2 is measured by a resolver rather 
than a rotary encoder as used on joint 1. This resolver directly measures the angular 
position and velocity of the circular spline with respect to ground through a gear pair with a 
6.67:1 reduction ratio. 
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Figure 2.5: Joint 2 sensor locations 
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The third and final testing station, as illustrated in figure 2.6, relies on a design 
configuration substantially different from the two other joints. In particular, while the 
wave-generator is still driven by a motor mounted to the circular spline, the output rotation 
of the joint is transmitted through the flexspline and the circular spline is attached to 
ground. This unique design not only allows for unrestrained rotation of the output link but 
also provides a facility for examining harmonic-drive behavior under a different operating 
configuration. From the sensor locations outlined in figure 2.7, it can be seen that, unlike 
the other two joints, the input encoder directly measures absolute motor rotation and there 
is no input tachometer. On the output side, a resolver measures the flexspline rotation with 
respect to ground through a gear-pair with a 7.2:1 reduction ratio. Since the torque sensor 
connects both the motor and circular spline to ground, it can measure directly the combined 

torque of the wave-generator and circular spline, or equivalently, the torque seen at the 
flexspline. 



Figure 2.6: Cross-section of joint 3 
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Figure 2.7: Joint 3 sensor locations 


2.1.2 Sensor Processing and Interpretation 

In order to monitor sensors and control motors on each joint, the data-acquisition 
system and motor-controller illustrated in figure 2.8 was used. Using a software 
environment called CONDOR, developed at the MIT Artificial Intelligence Laboratory, 
computer code was compiled on a Sun 3/180 workstation and loaded onto a series of 
Ironies processor boards. When executed, this code could then collect sensor information 
from the data-acquisition boards and sensor-processing circuits shown in figure 2.8. To 
control the motors, desired motor currents were specified on the D to A board and then 
amplified through an Aerotech 4020-LS DC Servo Amplifier. By monitoring the amplifier 

current with the A to D board, measurements of motor current and torque could be 
gathered. 

In order to extract useful information about harmonic-drive behavior, a thorough 
understanding of sensor measurement is required. First, relationships must be made 
between the measurements seen by the sensors on each joint and the rotations and torques 
on the individual harmonic-drive components. Given the variable names defined in table 
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2.3, the input and output rotations with respect to ground are defined in table 2.4, and the 
sensor-measurement relationships are defined in table 2.5. Second, as presented in table 
2.6, sensor resolution and accuracy must be characterized, and, lastly, sensors must be 
calibrated. Experimental calibration procedures and results are reported in Appendix A. 
From this catalog of sensor information, experimental data can now be used to accurately 
measure the dynamic, static, and kinematic performance of the three harmonic-drive test 
specimens. 


Table 2.3: 

Description of important harmonic-drive variables 

®wg, Gfs, and ^cs 

Angular position of the wave-generator, flexspline, and circular spline, 
respectively, as referenced from ground. 

^wg, ®fs, and ^cs 

Angular velocity of the wave-generator, flexspline, and circular spline, 
respectively, as referenced from ground. 

Twgi Tf s , and T^ 

Torque experienced by the wave-generator, flexspline, and circular 
spline, respectively. 


Table 2.4: Absolute joint rotation in terms of harmonic-drive variables 
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2.2 Kinematic Error Measurement 

The goal of my experimental investigation of harmonic-drive kinematic properties 
was to accurately characterize the position-error signatures of the three harmonic-drive 
specimens. As a result of this investigation, I not only confirmed the expected position- 
error magnitude and frequency distribution over typical frequency ranges but also 
characterized transmission accuracy at other rotational frequencies. Based on high- 
resolution experimental findings, evidence was found to trace the origins of kinematic 
inaccuracy to flexspline and circular-spline gear error. This section outlines the method 
used to collect kinematic-error data as well as the results and conclusions gained from this 
experimental analysis. 


2.2.1 Experimental Procedures 

Due to complications in separating harmonic-drive position-error components from 
kinematic inaccuracies inherent in the test apparatus, two different experimental procedures 
were required to measure position-error data. The general procedure was used for the 
majority of experimental trials while the phase-shift procedure was required to resolve 
transmission error that varied over one output rotation of the harmonic-drive. These two 
procedures are outlined below. 


2.2.1.1 General Procedure 


In order to collect position-error data, the three harmonic-drive testing stations were 
used to measure input and output rotation over a range of motion. By assuming quasi¬ 
static conditions, kinematic-error data was collected while the joint was rotated at a small 
but non-zero velocity. The typical procedure for data-acquisition and processing is listed 
below. 




By sending a constant current to the joint motor that exceeds the transmission 

starting torque, rotate the harmonic drive at the slowest velocity at which 
resonance vibration is minimal. 

Collect input and output position information from the joint sensors at equally 
spaced time intervals over a given motion range. 
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1 


0 
0 
0 

By adjusting the range of transmission rotation and the number of data points 
accumulated, FFT spectra spanning a wide range of frequencies were collected. Assuming 
that at least two data points per cycle are required to measure position-error frequency 
components, position-sensor resolution limits the cycles per angular position, or spatial 
frequency, that can be measured on each harmonic drive. These spatial frequency upper- 
bounds are listed in table 2.7. 




. 
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output encoder 

250.0 


output resolver 

1365.0 

IB11II!!!I!!!!11 

input encoder 

1000.0 


Calculate the position error, 0 err , using the equation: 


e er r = 


6 


in 


gear ratio 


-6 


out 


( 2 . 1 ) 


Plot the resulting position-error signal versus the number of input, or wave- 
generator, revolutions to analyze the shape of the error signature. 

Linearly interpolate the time-based position-error signal to produce data points 
that are equally spaced by input rotation. 

Take an FFT of the interpolated data vector to identify the frequency 
components of the error signature in terms of input revolutions. 


2.2.1.2 Phase-Shift Procedure 

Since a typical harmonic-drive position-error signature is composed of waveforms 
located at distinct spatial frequencies, most error components arising from harmonic-drive 
inaccuracies could be easily distinguished from inherent kinematic errors in the testing 
apparatus. Unfortunately, since position inaccuracy in the output sensor gear appeared at 
the same frequency as a potential harmonic-drive error component, a separate data- 
collection procedure was required to separate these two error components. More precisely, 
using the data-collection procedure above, harmonic-drive position-error components 
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occurring once every output revolution, and harmonics, were contaminated by inaccuracies 
in the output sensor gear which appeared at the same frequency. To distill the contribution 
of the harmonic drive from this aggregate error waveform, a second set of position-error 
data in which the harmonic drive was physically rotated by a fraction of a revolution 
relative to the output sensor gear was required. In other words, after the first set of 
kinematic-error data was obtained, the angular position of the harmonic-drive relative to the 
output gear was shifted by a known amount before the second data set was collected. By 
taking an FFT of these two signals, subtracting them, and scaling the result appropriately, 
the harmonic-drive error-component could be accurately identified. 

The theoretical basis for this experimental procedure can be clarified by assuming 
that the position-error waveform collected from the harmonic drive testing apparatus has 
three components: (1) h^G), a first error signal at a spatial frequency of once-per-output 
revolution , (2) h 2 ( 0 ), a second signal at a spatial frequency of once-per-output revolution, 
and (3) g(0), a third signal composed of various other components at different frequencies, 
where 0 is a measurement of angular position. Therefore, the aggregate position-error 


signal, h suml (0), can be defined as 

hsum i(e)=hi(e)+h 2 ( 0 ) + g( 0 ) ( 2 . 2 ) 

Now assume that a second signal is collected in which the first component was rotated by a 
known phase shift, A0, in the positive 0-direction. This signal, h sum2 (0), is described by 

hsum 2(0)=h!(0-A0)+h2(0) + g(0) (2.3) 

Taking the Fourier transform of both equations yields 

Hsumi(o)) = Hi(co) + H 2 (cd) + G(go) , and (2.4) 

Hsum 2 (co) =Hi(co)e- lCoA0 + H 2 (ct)) + G(go) (2.5) 


where CD is the spatial frequency in radians per rotational unit. By subtracting these two 
equations, it can be seem that the H 2 (cd) and G(co) terms disappear leaving an equation in 
terms of H^co) and the two aggregate error signals: 

Hi(CD) " H SU m2(C0) 

l _ e -ico A9 ’ { • ) 
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for co A0 * 2k n , n = 0, 1,2,... . Therefore, if a phase shift, can be accurately introduced 
into the harmonic-drive error component, it can be easily extracted from two sets of 
position-error data using this equation. 


The procedure used to collect and process two sets of position-error data according 
to this theory is itemized below. 



Set the zero setpoint for input and output position sensors on the given testing 
station. 


a 


Move the joint at a slow, non-resonating velocity by sending a current to the 
joint motor. 

Read input and output position data at equally spaced time intervals. 
Calculate the kinematic error from equation (2.1). 


@ Filter this position error data at every time step using a linear FIR digital low- 

pass filter to prevent aliasing of high-frequency position error components 
into sampled data. 

[6] Return the joint to its original position after this data run is complete. 


Interpolate the first set of position error data to generate data points that are 
equally spaced by input rotation rather than time. 


0 


Introduce a phase shift into either the harmonic drive, the output gear, or both 
by rotating that component a measured number of degrees. 

Collect the second set of data using the same procedure used to gather the first 
set. 



Interpolate this data set over the identical range of rotation used in the first set 
of data. 


Take an FFT of both data sets to yield frequency-domain data. 

Subtract the two resulting vectors of frequency domain data at all frequencies 

and scale the difference according to equation (2.6) to yield the magnitude 
spectrum of the phase-shifted position-error component. 

Take an inverse FFT of this spectrum to produce the waveform of the phase- 
shifted position-error component if desired. 
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2.2.2 Discussion of Results 

The results of several experimental tests on harmonic-drive kinematic error are 
presented in this section. After the importance of these results has been analyzed, 
comparisons will be made to catalog predictions of harmonic-drive position-error, and 
modeling representations of observed behavior will be discussed. Based on these findings, 
conclusions will be made about kinematic error in harmonic drives. 


2.2.2.1 Characteristic Harmonic-Drive Error 

For organizational clarity, I have divided my position-error experimental findings 
into four sub-sections. The first section describes the most significant harmonic-drive 
error, which appears at frequencies corresponding to multiples of input, or wave-generator, 
rotations. The second section presents results from high-resolution tests on these same 
error components. The third section focuses on error components that vary with output 
rotation while the fourth section focuses on high-frequency error signals which oscillate 
many times every wave-generator revolution. In each of these sections, experimental data 
will be interpreted and factors which can influence measured results will be discussed. 

2.2.2.1.1 Error at the Input Rotation Frequency 

Figures 2.9 through 2.14 illustrate the typical position-error waveforms observed 
on the three harmonic drives as well as the spatial FFT results which indicate the frequency 
content of these error signatures. For all of these plots, position-error magnitude is 
measured in degrees of output (flexspline or circular spline) rotation and plotted in terms of 
input (wave-generator) revolutions. In both the FFT and waveform plots, significant 
energy content can be seen at frequencies lower than one cycle-per-input-revolution. These 
low-frequency error components are due primarily to position inaccuracies in the output 
sensor gear train and should not be mistaken for harmonic-drive error. Unfortunately, the 
tooth-passing frequency of the output sensor gear-train appears on the joint 2 FFT 
spectrum at 3.0 cycles-per-input-revolution. When viewing this plot, the error component 
observed at this frequency should not be attributed to the harmonic-drive transmission. 

Table 2.8 summarizes the peak-to-peak error magnitudes observed on the three harmonic 
drives. 
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Table 2.8: Total peak-to-peak position error 

.. 

for the three harmonic drives 

mm M Wm iws »; :& sssg: . : 

v.v.v.v 

iii! 

XvXvX 



lIlBillH 


•x’:*:’:*X" 

Xvx*:*x 

0.045 

0.030 

0.060 

•x , x*x‘xtt , x , x , x*x , x*x : x : xsx%%¥!w;wjx:y;"ftx:»:x:x"t*x , x , x , x , x | x*x , x‘xx*x‘ 

x-:wx*X‘X*x-x-x*XixxX‘X*x*:’X-x*x*x-X‘X*x->x-: : :-: : x:#:X;x-:-x*:-XvXxx 

Ouaiact value pilpiilill 


0.042 

0.053 

0.095 


As past research suggests, the kinematic-error results display a primary error 
component at 2.0 cycles-per-input-revolution with smaller contributing components 
appearing in the next eight harmonics. As noted by Nye and Kraml, [24], in many FFT 
spectra, such as the ones for joints 2 and 3, the error component of the last observable 
harmonic, eighteen cycles-per-input-revolution, was surprisingly large in amplitude. 
Harmonic content contributed typically less than sixty percent of the total error signal. 
Additionally, in many FFT spectra, such as the ones for joint 2 and 3, error components 
were seen at 1.0 cycles-per-input-revolution and subsequent harmonics. This observation 
is not surprising since any misalignment or improper mounting of harmonic-drive 
components could likely produce kinematic error at this frequency. 

Through qualitative observations made during experimentation, I noticed a strong 
sensitivity of error magnitude and composition on joint velocity. More specifically, slight 
changes in the joint velocity during data collection caused noticeable variation in the 
position error signatures. For example, some tests, especially on joint 2 and 3, revealed 
discrepancies in error amplitude of up to 100 percent. Through investigations into the 
dynamic behavior which will be discussed later, it was discovered that, due to the high 
fluctuating torques excited by transmission error during rotation, significant torsional 
deflection can occur in the harmonic drive which increases or decreases position error. 
Since the harmonic drives on joint 2 and 3 have lower stiffness than the one on joint 1, it is 
likely that they are more susceptible to kinematic-error contamination by torsional 
deflection. The results presented above were collected carefully at slow velocities which 
exhibited minimal resonance vibration. Despite this experimental prudence, the true 
accuracy of the above results can only be determined by understanding the dynamic 
behavior of the transmission. Consequently, due to this breakdown of the quasi-static 
assumption, future position-error testing should be conducted on a stationary transmission 
in order to ensure measurement accuracy. 

By disassembling and reassembling the harmonic drive testing equipment between 
experiments, I observed some interesting behavior in the position-error profiles. In 
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particular, for a few tests, the largest position-error component appeared at four cycles-per- 
wave-generator revolution rather than two, and harmonic content of the aggregate error 
waveform increased substantially. Fortunately, no significant changes in the overall 
position-error magnitude were observed in the testing apparatus due to joint reassembly. 
However, as Nye and Kraml note in [24], since alignment can influence the magnitude of 
kinematic error in harmonic drives, I would not have been surprised if less careful 
reassembly resulted in different kinematic-error measurements. 

2.2.2.1.2 High-Resolution FFT Analysis of Primary Error Component 

To satisfy my curiosity, I performed a few tests in which I increased the FFT 
resolution around the two-cycles-per-input-revolution frequency by increasing the number 
of sampled data points. Not surprisingly, I discovered that separate and distinct error 
components appeared at the primary transmission-error frequency. Figure 2.15 shows an 
FFT magnitude spectrum for joint 1 magnified around the frequency of 2.0 cycles-per- 
input-revolution. As illustrated in this plot, an error component not only appears at the 
frequency of 2.0, but also at 2.0125 cycles-per-input-revolution. These two components 
can be explained by assuming that most harmonic-drive inaccuracy occurs due to gear-tooth 
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placement error and realizing that, during joint rotation, the flexspline moves relative to the 
circular spline. Specifically, since the input position sensor on joint 1 measures wave- 
generator rotation with respect to the circular spline, one-half revolution of the wave- 
generator will rotate the tooth-meshing zone exactly 180 degrees around the circular spline 
and complete one-cycle of the circular spline gear-error waveform. At the same time, since 
the circular spline will rotate one tooth ahead of the flexspline after this 180 degree wave- 
generator rotation, the flexspline gear-tooth error waveform will have completed one 
complete cycle plus one gear tooth or (1 + (1/N)) cycles. Therefore, since the gear ratio, 
N, on joint 1 is 160, for every complete rotation of the wave-generator the circular spline 
gear-tooth error will have completed two complete cycles and the flexspline gear error will 
have completed 2(1+ (1/160)) or 2.0125 cycles. From this explanation, it can be seen that 
the error component in figure 2.15 that appears at two cycles-per-input-revolution is due to 
the circular spline gear-tooth-placement error and the error component that appears at 
2.0125 cycles-per-input-revolution is due to the flexspline gear-tooth-placement error. 
Figure 2.15 also illustrates an error occurring at 2.00625 cycles-per-input-revolution which 
is probably due to the non-sinusoidal nature of the other two neighboring components. 

This high-resolution FFT analysis of the characteristic harmonic-drive kinematic- 
error frequencies may prove useful for many purposes. First, the presence of two distinct 
error components is strong evidence to reinforce the belief that the primary source of 
kinematic error in harmonic drives is due to gear-tooth inaccuracy in the flexspline and 
circular spline. Additionally, since this method can be used to extract quickly the individual 
position-error components due to the flexspline and circular spline, it may prove useful for 
characterizing the manufacturing accuracy of these two components as well as their 
individual performance accuracy under varied operating conditions. Lastly, since the 
amplitude modulation, or beating, that can appear in kinematic-error signatures is a product 
of these two distinct error components, the amount of beating in a given drive can be 
identified easily from the magnitude of these two individual errors. 


2.2.2.1.3 Error at the Output Rotation Frequency 

Using the phase-shift procedure described previously to separate the harmonic- 
drive error components at low spatial frequencies from the position inaccuracies of the 
output sensor gear-train, the FFT spectra shown in figures 2.16 and 2.17 were found. For 
both of these plots, position error was measured in degrees of output position with respect 
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Figure 2.16: FFT of joint 1 position error in cycles per output revolution 
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Figure 2.17: FFT of joint 3 position error in cycles per output revolution 
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to the frequency of output (circular-spline or flexspline), rather than input, rotation. From 
figure 2.16, it can be seen that the joint 1 harmonic drive exhibits a low-frequency error 
component occurring once per output revolution. The FFT peak that appears at 8.0 cycles- 
per-output-revolution is caused by inaccuracy in the output sensor pinion and should be 
ignored. In Figure 2.17, the error components for the joint 3 harmonic drive appear 
primarily at 1.0 and 3.0 cycles-per-output-revolution. Again, the FFT spike occurring at 
7.2 cycles-per-output-revolution is due to the output sensor pinion and should not be 
mistaken for harmonic-drive error. In both of these plots, the aggregate low-frequency 
error component due to harmonic-drive transmission-error remains significantly less than 
the harmonic-drive error observed at multiples of the input rotation frequency. 
Unfortunately, due to the limitations of the testing apparatus, the low-frequency position 
error of the joint 2 harmonic drive could not be measured. 

By performing a theoretical sensitivity analysis of the data-collection procedure, I 
determined that deviations in the final FFT error results due to known measurement 
uncertainly in the apparatus should be less than five percent. However, the results from 
several different trials showed variations in harmonic-drive error magnitude approaching 
sixty percent. This variance illustrated the high sensitivity of the low-frequency harmonic- 
drive error on operating conditions such as rotational velocity. 

The source of this low-frequency harmonic-drive error is probably a combination of 
harmonic-drive gear error and assembly errors. For example, it is conceivable that tooth- 
placement errors on the rotating output harmonic-drive component (circular spline or 
flexspline) could causes position deviations over one output rotation. Additionally, any 
misalignment of the flexspline or circular spline could cause similar position errors at 
multiples of the output rotation frequency. 


2.2.2.1.4 Error at High Frequencies 

Another set of experiments that I performed focused on higher rotational 
frequencies in order to detect the passing of harmonic-drive teeth in the position-error data. 
As figures 2.18 and 2.19 illustrate, essentially no position error components are visible at 
high spatial frequencies. Since the joint 2 harmonic drive has 320 gear teeth which engage 
for every input revolution, any contribution of the tooth-profile error to the aggregate 
position error should appear in figure 2.18 at a frequency of 320 cycles-per-input- 
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revolution. Similarly, since the joint 3 harmonic drive contains 160 gear teeth, the tooth¬ 
passing frequency should be seen at 160 cycles-per-input-revolution in figure 2.19. 
However, for both joints, no detectable error component was found at these frequencies. 
This unexpected result can be justified by recognizing that, at any point in time, the 
harmonic drive has several teeth in contact due to the concentric nature of its design and the 
shape of the wave-generator. Close observation revealed that the harmonic drives tested 
typically had up to twenty teeth in contact on both sides of the wave-generator. Given this 
result, it is unlikely that the individual error of a single gear tooth could manifest itself in 
the aggregate position-error spectrum. Unfortunately, the high-frequency error 
components on the joint 1 harmonic drive could not be measured due to resolution 
limitations of the output position sensor. 


2.2.2.2 Comparison to Catalog Values 

As listed in Table 2.8 in section 2.2.2.1.1, catalog predictions of harmonic-drive 
inaccuracy are compared to the combined amplitude of the error components measured at 
multiples of the input rotation frequency. From these values, it can be seen that the peak- 
to-peak error magnitude of the observed waveforms was less than or near the 
manufacturer’s quoted accuracy ratings. Also, as noted in the catalogs, the primary error 
component in all cases occurred twice every wave-generator rotation, but the catalogs failed 
to offer any predictions about the harmonic content of these signals or error components at 
other frequencies. 


2.2.2.3 Modeling Considerations 

Having measured the magnitude of the harmonic-drive kinematic-error components 
at a range of different frequencies, a mathematical representation of this behavior can be 
easily constructed. Specifically, by using a harmonic series, the aggregate effect of the 
important individual error components can be compiled: 

Gerfn = Ai sin (cOi0i + (jq) + A 2 sin (co 2 02 + <j>2) + .... (2.7) 

where A is the half-amplitude, co is the spatial frequency, 0 is a measure of rotation, and <J) 
is the phase of the given error component. Depending of the required accuracy of this 
representation, the number of terms in the harmonic series can be varied. Typically, in the 
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modeling efforts to be discussed later, I relied mainly on the two largest error components, 
which appeared at twice- and four-times every input revolution. The resulting harmonic- 
series representation of this error is 

Gerfn = Ai sin (29 in + <>!) + A 2 sin (40 in + (J> 2 ) , (2.8) 

where 0^ is the rotation of the harmonic-drive input. To be truly accurate, as indicated in 
the high-resolution tests presented above, two distinct error components should be included 
for each of the two terms in this series to account for the individual error contributions of 
the flexspline and circular spline which occur at slightly different frequencies. However, 
since it is unlikely that this small change will significantly influence performance of a 
dynamic model, this effect will be ignored. 


2.2.3 Kinematic Error Conclusions 

As discussed above, the experimental investigation of harmonic-drive kinematic 
error has confirmed previously observed behavior of position-error profiles as well as lent 
new insight about the sources and stability of harmonic-drive transmission inaccuracy. In 
particular, the magnitude of all observed position error-signals agreed with manufacturer’s 
specifications and the primary frequencies of aggregate error appeared at 2.0 cycles-per- 
input-revolution and subsequent harmonics. Additionally, these primary error components 
were confirmed to be produced by two separate error signals most likely due to the gear- 
tooth error on the circular spline and flexspline. Furthermore, a small but measurable 
amount of kinematic error was found at the frequency of output revolution while no 
significant position-error components were detected at the significantly higher tooth¬ 
passing frequencies. Based on these results, a simple harmonic series of the individual 
error components was generated in order to predict harmonic-drive kinematic error. Since 
transmission compliance can influence position-error measurements even at small rotational 
velocities, it is recommended that future kinematic tests on harmonic drives be performed 
under truly static circumstances. 
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2.3 Stiffness Testing 

The purpose of my investigation into harmonic-drive flexibility was two-fold: (1)1 
wanted to collect accurate stiffness data for the three harmonic-drive specimens to include 
in my dynamic models, and (2) I hoped gain a better understanding of possible ranges and 
sources of variation in typical stiffness profiles. As a result of my investigation, I not only 
characterized the harmonic-drive non-linear stiffness behavior and hysteresis loss, but also 
discovered that harmonic-drive stiffness measurements are surprisingly fickle and very 
dependent on assembly and operating conditions. This section is devoted to describing the 
testing procedure used to collect stiffness data as well as presenting the results of my 
investigation. 

2.3.1 Experimental Procedure 

The three harmonic-drive test stations were used to collect stiffness data by locking 
the wave-generator to the circular spline and applying a range of loads to a long output link. 
Figure 2.20 illustrates the testing equipment for joint 1 used to apply loads to the harmonic 


locked input 



56 
















Chapter 2: Stiffness Testing 


drive which is typical of the apparatus for the other two joints. The testing procedure used 
on each joint is as follows. 



Lock the wave-generator of the harmonic drive to be tested to its circular 
spline with a rigid fixture. 



Apply a load to the end of the output link of the joint using the screw-crank 
apparatus. 



Read and store output rotation data from the output position sensor and output 
torque data from the torque sensor. 



Continue collecting data for a range of positive and negative applied loads 
until the stiffness curve is complete. 


This straightforward procedure was used to collect multiple sets of stiffness data on each 
transmission and generate the results presented below. 


2.3.2 Discussion of Results 

The stiffness results collected on the three harmonic drives can be divided into three 
separate categories. First, a typical stiffness profile for each of the tested transmissions is 
discussed and variations in these measurements with environmental and operating 
conditions are outlined. Second, these stiffness results are compared to manufacturer’s 
quoted values in order to identify and explain discrepancies, and, lastly, theoretical 
methods for capturing the observed stiffness behavior are considered. 


2.3.2.1 Characteristic Stiffness Profiles 

Figures 2.21, 2.22. and 2.23 illustrate the measured stiffness curves on the three 
harmonic drives. Due to torque-sensor loading limitations, the joint 1 harmonic drive was 
loaded to 35 percent of its rated maximum torque, the joint 2 harmonic drive was loaded to 
45 percent of its rated maximum torque, and the joint 3 harmonic drive was loaded to 90 
percent of its rated maximum torque. These torque limitations, however, encompass the 
entire range of dynamic behavior seen by the harmonic drive in experiments to be discussed 
later. All three stiffness profiles exhibit noise due to limitations on position sensor 
resolution. 
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3 harmonic-drive stiffness curve 1 


As mentioned in manufacturer’s catalogs, all three harmonic drives exhibit a non¬ 
linear stiffness with appreciable hysteresis loss. As illustrated in figure 2.21, the harmonic 
drive on joint 1 displays a fairly well-behaved stiffness profile with a nearly linear stiffness 
and little hysteresis over the test range. As explained in the catalog, increased loading on 
this harmonic drive should increase the non-linearity of the stiffness curve. The stiffness 
profile of the joint 2 harmonic drive, as illustrated in figure 2.22, shows significant non¬ 
linearity despite being loaded over only 45 percent of its rated torque range. Additionally, a 
hysteresis loss is visible which is substantially larger than the loss observed on the other 
two harmonic drives. Since no backlash was present in the joint 2 testing apparatus, this 
hysteresis can be blamed on a combination of static friction and very low stiffness levels at 
low applied torques, also called soft wind-up, which can sometimes act like backlash. The 
stiffness curve for the joint 3 harmonic drive, as presented in figure 2.23, shows very 
typical behavior with a non-constant slope and a modest amount of hysteresis. 

Another set of stiffness tests that I performed were aimed at gaining a qualitative 
understanding of how stiffness varied with environmental and operating conditions. In 
particular, I used the joint 1 testing station to collect data for a range of different preloads 
and wave-generator angles. Figure 2.24 shows three stiffness curves generated from data 
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Figure 2.24: Joint 1 stiffness curves for different wave-generator angles 
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Figure 2.25: Joint 1 stiffness curves for different wave-generator preloads 
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collected with the wave-generator major-axis locked in different positions relative to the 
circular spline. As the difference between these curves indicates, variations in stiffness of 
up to 25 percent may be measured due to different orientations of the wave-generator. 
Figure 2.25 illustrates three different stiffness curves collected when the preload between 
the harmonic drive gear teeth was varied by changing the depth of the wave-generator in the 
flexspline. While wave-generator depth is not a direct indicator of gear-tooth preload, it 
was observed to influence the radial loading of the harmonic drive in many cases. The test 
labeled “high preload” corresponds to the case where the wave-generator was pushed into 
the flexspline to its maximum depth while the low-preload case saw the wave-generator at 
minimum depth. From the stiffness curves in this plot, it can be seen that slight variations 
in preload can also influence the value of harmonic-drive stiffness. 

In general, from the range of behavior illustrated in the stiffness results, it can be 
concluded that the variation in stiffness profiles between different harmonic drives can be 
substantial. As seen on joints 1 and 2, one drive can behave relatively linearly in the low- 
to-medium torque range while another can exhibit a strikingly non-linear profile over the 
same loading range. Additionally, the presence of hysteresis due to phenomenon such as 
friction and soft wind-up may be dominant in one harmonic drive while unnoticeable in 
another. Furthermore, the substantial variations in harmonic drive flexibility due to such 
factors as preloading and input orientation indicate the sensitivity of stiffness 
measurements. This gloomy forecast for the predictability of harmonic-drive flexibility 
grows more overcast when these results are compared to manufacturer’s quoted values. 


2.3.2.2 Comparison to Catalog Values 

In manufacturer’s catalog [4], a piecewise-linear-approximation method is proposed 
for predicting harmonic-drive stiffness. In particular, two-piece linear fits are quoted for 
the case of a standard stiffness profile and an optimal, or maximum-stiffness, profile. 
Figures 2.26, 2.27, and 2.28 show these two piecewise linear approximations compared to 
the measured stiffness for each harmonic-drive model. As these results illustrate, the joint 
1 harmonic-drive stiffness is considerably different from the quoted values while the 
transmission stiffness on joints 2 and 3 more closely resembles the catalog stiffness 
profiles. More specifically, on joint 1, the measured harmonic-drive stiffness is about 30 
percent greater than the optimal stiffness curve and about 60 percent greater than the 
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Figure 2.26: Joint 1 linear approximations to stiffness data 
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Figure 2.27: Joint 2 linear approximations to stiffness data 
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Figure 2.28: Joint 3 linear approximations to stiffness data 


standard stiffness curve. As noted in past research, this unexpectedly high stiffness can 
probably be explained by a high gear-tooth preload in the harmonic drive. The lack of 
hysteresis loss observed on joint 1 also supports this postulation since a higher preload 
would eliminate any soft wind-up in the transmission. On joint 2, the measured stiffness 
while loading the harmonic drive is slightly lower than the standard stiffness approximation 
while the optimal quoted stiffness agrees closely with the unloading portion of the 
measured stiffness curve. On joint 3, both the loading and unloading portions of the 
measured stiffness curve have a slope that agrees with the standard stiffness approximation 
within about 10 percent. In partial confirmation of the discrepancies between predicted and 
measured stiffness values, manufacturer’s catalog [3] warns that a 30 percent variation 
between quoted and actual stiffness values is not uncommon. In terms of hysteresis, 
catalog [3] states that hysteresis loss rarely exceeds 0.033 degrees of output torsion and 
that soft wind-up is kept below 0.05 output degrees of rotation. This prediction holds true 
for joints 1 and 3 but is greatly underestimated for the joint 2 harmonic drive. 

This comparison between experimental data and catalog predictions exposes the 
shortcomings of manufacturer’s estimates. Multiple stiffness curves collected on each 
harmonic-drive assembly illustrated that experimental results were fairly repeatable to 
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within about ten or fifteen percent. However, measured values of stiffness and hysteresis 
loss differed unpredictably and significantly from quoted values. In a few instances such 
as the joint 1 stiffness profile and the joint 2 hysteresis loss, measured values even fell 
outside of the manufacturer’s predictability tolerances. 


2.3.2.3 Modeling Considerations 

Figures 2.29, 2.30, and 2.31 show cubic approximations to the stiffness curves for 
the three harmonic drives. In all cases, the cubic polynomial was fit to the portion of the 
stiffness data collected during joint loading rather than unloading. From these results, it 
can be seen that cubic approximations can accurately capture the shape of experimental 
stiffness data. Based on this result, cubic representations of non-linear harmonic-drive 
stiffness will be used in dynamic modeling efforts to be discussed later. 

One measure of the accuracy of harmonic-drive stiffness representation is whether 
or not it can replicate the observed hysteresis losses. Capturing these effects is dependent 
on type of dynamic model used for the harmonic-drive and the placement of the flexibility 
within this model. Conversely, given a harmonic-drive model which can describe 



Figure 2.29: Cubic-tit ot joint 1 harmonic-drive stittness data 
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Figure 2.30: Cubic-fit of joint 2 harmonic-drive stiffness data 
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Figure 2.31: Cubic-fit of joint 3 harmonic-drive stiffness data 
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hysteresis loss, the measured stiffness profiles can be used to make estimates of the friction 
present in harmonic drives. As discussed later in this document, dynamic models of 
harmonic drives were developed which capitalized on the friction information stored in the 
hysteresis loss of the stiffness measurements. 

As the observed variation in measured stiffness values indicates, a truly accurate 
model of harmonic-drive stiffness should describe the dependence of the stiffness profile 
on parameters such as input angle and gear-tooth preloading. Although I did not perform 
enough tests to develop a quantitative characterization of these effects, additional 
measurement of harmonic-drive stiffness over different operating conditions could be 
performed easily to characterize empirically an accurate stiffness function. Alternatively, a 
much more difficult but meaningful approach to this problem would be to develop an 
analytical model which can predict stiffness variations based on harmonic-drive geometry 
and the physics of operation. I will save this formidable task for later research. 


2.3.2.4 A Note on Input Diametral Stiffness 

If the wave-generator of an assembled harmonic drive is rotated manually, the 
torque required to rotate the transmission varies noticeably over one wave-generator 
revolution. If the circular spline is then removed and the same test is performed on the 
wave-generator and flexspline alone, the same torque variation is still observed. Although 
difficult to notice on some drives due to excessive friction, this cyclic torque has a tendency 
to rotate the wave-generator to preferred orientation when the load is released from the 
input. These qualitative observations indicate the influence of a diametral flexibility which 
acts on the harmonic-drive input due to manufacturing and assembly errors in the wave- 
generator and flexspline. Quick estimates of this input torque amplitude placed its 
magnitude at about 0.01 N-m for joint 1 and about 0.001 N-m for joints 2 and 3. For most 
applications, this cyclic torque is significantly less than the operating torque of the 
harmonic drive and can be ignored. However, in some applications, especially at low 
velocities, the presence of this cyclic torque can produce variations in rotational velocity at a 
frequency of twice-per-input-revolution. Consequently, to improve performance of a 
dynamic model, a simple sinusoidal torque function which varies at the this frequency can 
be included to replicate the effect of this input diametral flexibility. 
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2.3.3 Stiffness Conclusions 

The most significant conclusion that can be derived from the results in this section 
in that accurate measurement and prediction of harmonic-drive stiffness profiles can be very 
difficult. Not only are catalog estimates likely to be significantly different from actual 
values, but experimental measurement can also vary substantially as a function of 
operational conditions such as wave-generator angle and preload. To foreshadow issues 
which appear later in this document during dynamic modeling discussions, the static and 
coulomb friction which are manifested as hysteresis loss in stiffness measurements 
compound the problem of extracting accurate stiffness information. Since these friction 
phenomena are unpredictable and difficult to separate from pure transmission stiffness, the 
accuracy and usefulness of harmonic drive-stiffness measurements becomes even more 
nebulous. 


67 


CJia2t_er_Jj_Sj_arting-Torque Measurem ent 


2.4 Starting-Torque Measurement 

Due to static friction in the harmonic-drive transmission, a finite input torque is 
required to effect output rotation. The purpose of my investigation into this harmonic drive 
property was to obtain accurate values for the amount of static friction present in harmonic- 
drive transmissions and, more importantly, gain insight about the mechanism and behavior 
of velocity-independent friction in harmonic drives. As a result of this experimental 
investigation, predicted starting-torque magnitudes were confirmed but their substantial 
variation with several operating and environmental conditions exposed the complexity of 
physical mechanism behind the values. This section outlines the procedure used for these 
tests as well as the experimental observations and resulting conclusions. 


2.4.1 Experimental Procedure 

The procedure used for these experiments was designed to collect accurate 
measurements of starting torque in order to understand the amount of static friction present 
in a stationary harmonic drive. This procedure is listed below. 

Warm-up the given testing station to a typical operating temperature by 
running the harmonic drive for several minutes. 

Gradually apply a current to the input motor while collecting input torque data 
from the current sensor, and position data from the input encoder. 

Calculate the harmonic-drive starting torque by subtracting the catalog value 

for the motor starting torque from the maximum motor torque reached just 
before the joint began to rotate. 

starting torque in each harmonic-drive testing station is composed of the 
individual starting torques due to (1) the input motor, (2) the harmonic drive, and (3) the 
output bearings. Since the starting torque exerted by the output bearings, as measured 
from the input, will be very small, accurate values for the harmonic-drive starting torque 
can be found simply by subtracting the motor starting-torque from the measured value. 


0 

0 

0 

The total 


68 



C]uigie^2j_Jiia^ M easurem ent 


2.4.2 Discussion of Resuits 

The results of several qualitative and quantitative investigations into starting-torque 
behavior are presented below. First, typical starting-torque values will be discussed in 
light of environmental and operational conditions which can influence their magnitude. 
Additionally, these values will be compared to catalog predictions and analyzed to gain 
insight about potential modeling methods. Based on these results, conclusions about 
harmonic-drive starting-torque behavior will be drawn. 


2.4.2.1 Characteristic Starting-Torque Behavior 

Table 2.9 lists the measured starting torques for all three harmonic drives. These 
starting-torque measurements, however, were observed to vary considerably with the 
rotational position of the harmonic-drive output (flexspline or circular spline). These 
variations, as listed in table 2.9, can be explained by noting that any misalignment of the 
flexspline or circular spline, for example, is likely to influence gear-tooth meshing and, 
consequently, gear-tooth friction over one output revolution. Since these starting-torque 
values provide a direct measurement of static friction, the substantial variation in starting 
torque with output position illustrates the sensitivity of static friction in harmonic drives. 
To further support this unfortunate consequence, additional observations revealed a strong 
dependency between harmonic-drive starting torque and operating temperature. In a few 
cases, if testing equipment was not warmed-up before data collection, starting torques 
could increased by up to 80 percent. This troublesome effect can be explained partially by 
the increase in lubricant viscosity which increases interface friction at lower temperatures. 
Additional variables which can constrain the accurate measurement and prediction of 
harmonic-drive static friction are wave-generator orientation and gear-tooth preload. As 
noted in previous research, since the orientation of the wave-generator as well as the 
applied preload can influence the tightness-of-fit between the flexspline and circular spline, 
static friction has been observed to waver considerably under the influence of these two 
factors. For example, experimental tests presented in [65] showed that starting torque 
could vary as much as 50 percent over one rotation of the wave-generator. Due to the 
limitations of the testing equipment, I was not able to confirm this result, however 
qualitative observations clearly demonstrated the importance of input angle and preload on 
starting-torque measurement. 
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Table 2.9: Input starting-torque values for the three harmonic drives 
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Measured value (N-m) 

0.084 

0.031 

0.005 

MgMIgillMlgll 

0.078 

0.035 

0.007 

Variation with output angle 

10% 

25% 

15% 


2.4.2.2 Comparison to Catalog Values 


In addition to listing the measured harmonic-drive starting-torques, table 2.9 


presents predicted values as quoted in manufacturer’s catalog [4]. By comparing these 
results, it can be seen that the predicted no-load starting torques agree fairly well with the 


experimental results. However, as discussed above, the accuracy of both the experimental 
and catalog values is questionable due to the large number of factors which can greatly 


influence the magnitude of harmonic-drive starting torque. 


2.4.2.3 Modeling Considerations 

Since, by definition, static friction occurs only when velocity is zero, developing an 
accurate model of this phenomenon is not critical for predicting harmonic-drive behavior at 
non-zero velocities. However, by assuming that velocity-independent friction, such as 
Coulomb friction, behaves similarly to static friction, conclusions can be drawn about 
friction components at non-zero velocities. For example, due to the preload between the 
harmonic-drive gear teeth, a non-zero contact and rubbing force will always exist. From 
the definition of Coulomb friction, this non-zero normal force should always generate 
proportional frictional losses at all speeds of harmonic-drive operation. A truly accurate 
model of this coulomb friction behavior should use the harmonic-drive geometry to derive a 
value for the gear-tooth preload which can then be used to calculate the resulting Coulomb 
friction. Alternatively, for the sake of simplicity, these frictional losses due to preloading 
can be approximated by either a constant frictional torque or a function which replicates the 
typical behavior of the Coulomb friction. By assuming that the static friction results have a 
magnitude and behavior similar to this velocity-independent friction function, estimates of 
harmonic-drive frictional losses at non-zero operating velocities due to Coulomb friction 
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can be made. These insights about velocity-independent friction gained from starting- 
torque measurements will also be supplemented by the dynamic damping measurements 
presented below. 


2.4.3 Starting-Torque Conclusions 

As illustrated by experimental results, harmonic-drive no-load starting torques can 
be relatively large in magnitude and very dependent on operating and environmental 
conditions. While experimental measurements agreed fairly well with catalog values, the 
dependence of these values on factors such as gear-tooth preload, input and output 
orientation, and operating temperature undercuts the validity and usefulness of these 
numbers. However, by assuming similarities between this observed static-friction 
behavior and inherent Coulomb friction in the transmission, a useful understanding of 
velocity-independent friction behavior at non-zero velocities can be gained. 


71 


Chapter 2 : Dynamic Friction Measurement 


2.5 Dynamic Friction Measurement 

All harmonic drives exhibit a finite power loss during operation due to gear-tooth- 
rubbing friction and damping in the wave-generator bearing. The purpose of my 
investigation into these harmonic-drive frictional effects was to obtain accurate 
measurements of the losses in the three harmonic-drive specimens as well as gain a better 
understanding of the origins and magnitudes of harmonic-drive friction. As a result of my 
experimental investigation, I not only estimated values for the velocity-independent offset- 
torque experienced by harmonic drives at non-zero velocities, but I also characterized the 
non-linear damping profile for each harmonic-drive specimen. Additionally, due to torque 
fluctuations excited by system resonance, I discovered that frictional losses in certain 
operating regimes can increase dramatically. This section is devoted to outlining the 
procedure used for these tests as well as presenting interesting experimental results. Based 
on these results, conclusions are drawn about the nature of frictional losses in harmonic- 
drive transmissions. 


2.5.1 Experimental Procedure 

The procedure used for these experiments was designed to collect dynamic friction 
data for a range of rotational velocities. This procedure is listed below. 

Warm-up the given testing station to a typical operating temperature by 
running the harmonic drive for several minutes. 

Send a constant current to the DC motor and wait until an equilibrium 
rotational velocity is reached. 

Measure the input motor torque at this velocity. 

Collect data at several different operating velocities by repeating the two 
previous steps. 

Assuming that the rotational velocity during data collection was truly constant, the input 
motor torque should exactly balance the damping torque at the given velocity. Therefore, 
by collecting input torque data over a range of different velocities, the damping profile for 
each harmonic drive was constructed. 


0 

0 

0 

0 
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2.5.2 Discussion of Resuits 

The results of my experimental investigation of harmonic drive dynamic losses are 
presented in this section. First, characteristic loss profiles will be presented and analyzed 
followed by a comparison of these values to catalog predictions. Based on this 
information, representations will be proposed to model friction behavior and conclusions 
will be made. 


2.5.2.1 Characteristic Friction Behavior 

By using the procedure described above, the input motor torque required to operate 
the harmonic drive at a constant velocity was used to determine the damping torque of the 
transmission at different velocities. Due to the harmonic-drive kinematic error, the 
equilibrium velocity reached by the harmonic drive in many trials experienced noticeable 
fluctuation. In these cases, the average equilibrium velocity was used to determine the 
damping torque. The resulting torque and velocity data that was collected on each testing 
station directly characterized the total friction torque in the system. This total damping 
torque included not only the friction present in the harmonic drive but also the damping 
present in the input motor and output bearings as well. By assuming that, since the 
harmonic drive output rotates at relatively small velocities, the damping contribution from 
the output bearings was negligible, the frictional torque due solely to the harmonic-drive 
transmission was calculated by subtracting the motor damping from the measured aggregate 
friction-torque. Figures 2.32, 2.34, and 2.36 illustrate the total measured friction-torque 
plotted with the motor damping, while figures 2.33, 2.35, and 2.37 show the extracted 
harmonic drive friction-torque component. 

From the six plots shown below, several observations about dynamic friction 
losses in harmonic drives can be made. First, as figures 2.32, 2.34, and 2.36 illustrate, 
the friction torque generated in the harmonic drive is much larger than any other damping 
losses in the testing equipment. This anticipated result pronounces the importance of 
accounting for transmission losses in harmonic-drive models and calculations. Second, as 
labeled in figures 2.33, 2.35, and 2.37, each drive passes through operating regions in 
which frictional losses are amplified due to dynamic resonance vibration. More 
specifically, at operating velocities where the position error of the harmonic drive excites 
joint resonances, the resulting torque fluctuations can increase friction in the transmission 
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Figure 2.37: Damping curve for the joint 3 harmonic drive 
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and dissipate extra energy. Depending on the location of these resonances as dictated by 
the design of the harmonic-drive system, the frictional losses excited in these regimes are 
likely to dominate harmonic-drive behavior for a significant range of operating velocities. 
Third, the shape of the damping curve and the torque-offset experienced at low rotational 
velocities can lend insight into the nature of the non-zero, velocity-independent friction- 
torque resulting from the harmonic-drive gear-tooth preload. As discussed in the previous 
section, Coulomb friction at the gear-tooth interface due to non-zero gear-tooth loading 
influences the starting-torque of the harmonic drive as well as friction-torque offset 
experienced at non-zero velocities. Rough estimates of this velocity-independent offset- 
torque are derived from the damping-curve profiles in figures 2.33, 2.35, and 2.37 and are 
labeled “constant friction torque”. As explained previously, these frictional offsets may be 
constant with respect to velocity but are likely to be influenced by other parameters such as, 
input and output rotational angle, operating temperature, and preload. 

As observed in most other harmonic-drive properties, dynamic damping is not 
immune to outside influence. Specifically, qualitative observations showed that damping 
could dramatically rise with an increase in gear-tooth preload. This unsurprising effect can 
undoubtedly be explained by higher gear-tooth rubbing losses induced by increased 
loading. As indicated in manufacturer’s catalogs, frictional losses typically increase with 
higher gear-ratios and lower operating torques and temperatures. Also, the type and 
amount of lubricant used in the transmission can strongly influence dynamic friction 
measurements. Lastly, although difficult to distinguish from velocity-independent 
frictional losses, it would not surprise me if dynamic losses experienced noticeable 
variation with output and input position. 


2.5.2.2 Comparison to Catalog Values 

Harmonic drive catalogs, [3, 4], characterize the dynamic losses experienced in the 
transmission by providing efficiency ratings at different operating velocities. These 
efficiency ratings are derived by comparing the friction torque in the transmission to the 
rated torque of the harmonic drive at different operating velocities. Harmonic-drive 
catalogs specify the rated torque for each drive by using bearing-life equations to calculate 
the operating torques required to ensure a constant wave-generator bearing life at all 
operating velocities. Using this definition of the rated torque and the experimental frictional 
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Figure 2.38: Joint 1 harmonic drive efficiency 



Figure 2.39: Joint 2 harmonic drive efficiency 
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Figure 2.40: Joint 3 harmonic drive efficiency 


losses, the resulting transmission efficiency was calculated for each harmonic drive. 
Figures 2.38, 2.39, and 2.40 illustrate these experimental efficiencies as well as the catalog 
efficiency curves from [4]. Given that the catalog values are provided as rough estimates, 
these three plots demonstrate that the catalog efficiencies are relatively similar to the 
experimental results. For all three experimental curves, the dramatic efficiency loss that 
occurs during resonance is apparent, and, since the location of these system resonances is 
dependent on harmonic-drive compliance as well as the inertias of the specific harmonic- 
drive system, it is not expected that catalog efficiency curves should capture this behavior. 
Additional discrepancies between experimental and predicted results can be accounted for 
by a number of environmental factors such lubrication or operating temperature. 


2.5.2.3 Modeling Considerations 

In order to accurately model harmonic-drive friction at non-zero velocities, three 
separate friction components must be represented: (1) velocity-independent friction, (2) 
velocity-dependent friction, and (3) friction at resonance. The first of these three 
components is the non-zero Coulomb-like friction that occurs in harmonic drives due to 
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gear-tooth preload. In its simplest form, this component can be represented by a constant 
torque at all velocities, but, as discussed above, a truly accurate model should incorporate 
influences from such factors as input and output orientation and gear-tooth preload. The 
second friction component, namely non-linear velocity-dependent damping, can be better 
understood by ignoring the resonance losses in the damping profiles of figures 2.33, 2.35, 
and 2.37. The resulting damping curves, as shown in figures 2.41, 2.42, and 2.43, can be 
approximated accurately by cubic functions. Given these two models of harmonic-drive 
friction components, the representation will be complete if the behavior of resonance losses 
can be captured. Unfortunately, this last friction component is somewhat more difficult to 
model since it depends on the dynamic torques experienced by the harmonic drive. 
Consequently, a representation which produces these effects must be developed in 
conjunction with a complete dynamic model of the transmission as will be discussed later. 


2.5.3 Friction and Damping Conclusions 

From the results presented in this section, I can conclude that frictional losses in 
harmonic drives at non-zero velocities are characterized by (1) a velocity-independent offset 



Figure 2.41: Cubic fit of joint 1 non-resonance harmonic-drive damping 
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torque, (2) a non-linear damping torque, and (3) an additional torque loss during system 
resonance. Useful measurements for the first two of these frictional components can be 
made by analyzing the shape of the damping profile of a given harmonic drive while 
accurate representation of the third friction mechanism is more elusive . Due to the 
dependence of transmission friction on environmental and operating conditions, catalog 
efficiency curves can only be used to get a rough approximation of harmonic-drive losses at 
different velocities. 
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2.6 Dynamic Response Measurement 

With the insight gained from the kinematic and static experimental results presented 
above, the dynamic behavior of harmonic-drives can now be analyzed and interpreted. The 
purpose of this final investigation is threefold: (1) to fortify an empirical understanding of 
harmonic-drive dynamic behavior by identifying relationships between the observed 
dynamic response and known transmission properties such as stiffness, kinematic error, 
and friction, (2) to generalize about typical harmonic-drive dynamics based on experimental 
information collected from the three harmonic-drive specimens, and (3) to fully characterize 
the dynamic step-response behavior of the three harmonic drives in order to establish a 
reference point for the generation of theoretical models. The results of this investigation 
illustrated that there are several dynamic properties common to all of the harmonic drives 
tested. In particular, all transmissions displayed dramatic non-linear resonance behavior 
when the spatial frequency of the kinematic error harmonics coincided with the first 
resonant mode of the harmonic-drive system. In all cases, these resonances appeared over 
a substantial portion of the harmonic-drive operating range and behavior in these areas was 
uncommonly turbulent and unpredictable. This section is devoted to describing the 
experimental procedure used on the three testing stations as well as analyzing the collected 
dynamic-response data. Based on the information interpreted from the dynamic-response 
plots, modeling issues will be addressed and conclusions will be drawn. 


2.6.1 Experimental Procedure 

In order to investigate the dynamic behavior of harmonic drives, step-response data 
was featured exclusively as an indicator of overall harmonic-drive performance. 
Specifically, data was collected on all joints for harmonic-drive rotation resulting from a 
series of step-commands in current sent to the input motors. As illustrated in figure 2.44, 
each testing station was configured with a different output inertia representing a typical load 
for the given harmonic-drive. Additionally, the input inertias varied as dictated by the 
motor armature and wave-generator on each joint. A listing of all these inertia values is 
given in table 2.10. Using these three testing facilities, the experimental procedure listed 
below was performed. 

Q] Warm-up the given testing station to a typical operating temperature by 
running the harmonic drive for several minutes. 
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Position the joint output at a known starting position. 

Send a step-command of current to the DC motor amplifiers. 


Collect data from the position sensors, current sensor, and torque sensor for 
several seconds of joint operation. 

Derive velocity and kinematic error response from the position-sensor data. 

Collect data at several different input-current step-commands by repeating this 
procedure. 


The motor-current commands used in each trial ranged from the lowest current that could 
overcome transmission starting-torque to the current-saturation limit of the motor 
amplifiers. The resulting velocity effected by these current commands spanned the typical 
recommended operating range for the three harmonic drives. In all cases, current step- 
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Table 2.10: Inertia parameters for th€ 

dynamic-response tests 





llliilliill inertia 

1.8e-4 kg-m 2 

4.8e-5 kg-m 2 

3.2e-6 kg-m 2 


3.8e-5 kg-m 2 

2.8e-5 kg-m 2 

2.0e-6 kg-m 2 


total Input inertia 

2.2e-4 kg-m 2 

7.6e-5 kg-m 2 

5.2e-6 kg-m 2 

output mass 

2.828 kg 

0.7027 kg 

0.7027 kg 

mounting arm longtn H** 

0.308 m 

0.282 m 

0.112 m 

total output Inertia 

0.268 kg-m 2 

0.056 kg-m 2 

0.0088 kg-m 2 


commands were incremented at regular intervals small enough to capture the variations in 
dynamic behavior with increased current. From the resulting myriad of data sets, a handful 
were extracted which effectively characterized the dynamic response over the entire 
experimental range of operation. 


2.6.2 Discussion of Results 

In this section, several time-response plots will nucleate discussion of harmonic- 
drive dynamic behavior. From these plots, several characteristic behavioral patterns will be 
highlighted and used to gain insight into possible harmonic-drive modeling representations. 
Lastly, based on this improved understanding of harmonic-drive dynamic behavior, 
conclusions will be made about typical and predictable harmonic-drive dynamic response. 


2.6.2.1 Characteristic Dynamic Behavior 

All of the results from my experimental investigation of harmonic-drive dynamic 
behavior are illustrated in a series of time-response plots featured in Appendix D. These 
plots illustrate the position, velocity, and torque information gathered for the three 
harmonic drives for several different step-responses which fully capture the range of 
observed transmission behavior. For convenience, table 2.11 lists the figure numbers for 
all of the experimental plots in Appendix D. Throughout this section, references will be 
made to several of these plots in order to outline key features of the observed dynamic 
response. The discussion of characteristic harmonic-drive dynamic response will be 
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Table 2.11: Summary of 

experimental data-plot numbers in 

Appendix D 






PXw*SMfiKi»Jw25;7K"R5>>>w*K*%v»;»v«v.v.y 

❖X*X‘X*X‘X*X*X*:‘Xy:XxX:X:X:X:X!X£X?XrX?X:X 

miiiiliiiii 

dolnt 3 


D.1.1 

D.2.1 

D.3.1 


D.1.2 

D.2.2 

D.3.2 


D.1.3 

D.2.3 

D.3.3 

Output Position 

D.1.4 

D.2.4 

D.3.4 

Output Velocity - view i 

D.1.5 

D.2.5 

D.3.5 

•x<^x*x , >x< , >>>x , >X , X , x , x , >>X*>x*>!*x*£S£%x*» , x*x , x , x*x , x , x»x , x , x , x , x , x , x*x*x*Xv 

Output Velocity 1 view 2 

D.1.6 

D.2.6 

D.3.6 

— 11MM—1 

D.1.7 

D.2.7 

D.3.7 

Dynamic Position error 

D.1.8 

D.2.8 

D.3.8 


structured by first outlining typical behavioral patterns as observed in the velocity response, 
and then discussing the influence of these patterns on torque, current, position and position 
error measurements. 


2.6.2.1.1 Typical Dynamic Behavior Patterns 


In all of the step-response plots collected, underlying commonalities in dynamic 
behavior can be seen. However, the most illustrative information about the nuances of 
harmonic-drive dynamic behavior can be found in the velocity time-response plots in 
Appendix D. Specifically, from the input velocity curves illustrated in figures D.1.2, 
D.2.2, and D.3.2, four unexpected dynamic characteristics observed in all three 
transmissions become apparent: 


A substantial portion of the operating range of each harmonic drive is 
contaminated by resonance vibration. 





These regions of resonance act as energy sinks which constrain the increase 
of average rotational velocity for increasingly larger current commands. 

An unpredictable jump in velocity can frequently occur when the transmission 
harnesses just enough energy to push through a system resonance. 

Due to variations in frictional resistance over one output revolution of the 
transmission, significant variation in rotational velocity can result. 
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These surprising features which are common to the measured response of all three 
harmonic-drives will be discussed in order. 

From the listed time-response illustrations, it can be seen that the input velocity of 
the three harmonic drives passes through regions where significant resonance vibration can 
occur. These resonant regions occur when the frequency of kinematic-error fluctuations, 
which occurs at multiples of input, or wave-generator, rotational speeds, coincides with the 
primary temporal resonance of the harmonic-drive system. Since the harmonic-drive 
transmission flexibility behaves non-linearly, there is no distinct frequency of system 
resonance. Instead, as the torque-fluctuation amplitude increases, the effective stiffness of 
the transmission increases which results in an increase of the frequency of system 
resonance. This variation in temporal resonant behavior is manifested in a range of 
operating velocities which can excite system resonance. On joint 1, two separate 
resonances can be observed at input speeds of around 4000 and 8000 degrees/second. The 
vibration in the first of these two regions can be attributed to the kinematic-error excitation 
which occurs four times every input revolution and the second region can be blamed on the 
twice-per-input-revolution position-error component. As implied, the resulting frequency 
of velocity fluctuation that occurs at each of these resonances is identical to the rotational 
speed of the harmonic-drive input times the spatial frequency of the responsible kinematic- 
error excitation. On joint 2, the amplitude fluctuations at resonance are smaller than 
observed on joint 1, but regions of resonant behavior can still be found at about 10000 
degree s/second and , to a lesser extent, 20000 degrees/second. These two regimes can be 
attributed to the kinematic-error fluctuations occurring at 2.0 and 1.0 cycles-per-input- 
revolution, respectively. On joint 3, three resonant areas can be observed at the 
approximate input velocities of 5000, 10000, and 20000 degrees/second, and can be 
directly associated with the kinematic error occurring at 4.0, 2.0, and 1.0 cycles-per-input- 
revolution, respectively. The lack of any significant system resonance on joint 1 
corresponding to a frequency of once-per-input-revolution is not surprising since, as 
illustrated previously in figure 2.10, no appreciable kinematic error component was 
measured at that frequency. Additionally, kinematic-error harmonics appearing at 
frequencies higher than 4.0 cycles-per-input-revolution have no effect on any of the 
observed dynamic responses since the resulting excitation frequencies exceeded the system 
resonant regions at the slowest experimental operating velocities. 

An unfortunate consequence of the resonance experienced by the harmonic drives is 
the dramatic dissipation of energy which occurs in these regions. Specifically, from figure 
D.1.2, it can be observed that input current commands from 2.0 to about 2.4 amps achieve 
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roughly the same final rotational velocity at the first region of resonance, and current 
commands ranging from 2.8 to 4.0 amps all produce the same final velocity in the second 
resonance regime. Similarly, as illustrated in figures D.2.2 and D.3.2, identical behavior 
can be observed around the resonance regimes on joints 2 and 3. The velocity limitation at 
resonance can be attributed to the dramatic increase in Coulomb-like friction at system 
resonance due to the large torque fluctuations experienced in the transmission. In other 
words, as a step command in current brings the system closer to resonant velocity, the 
torque fluctuations due to resonance begin to increase. With this increase in torque 
amplitude, energy is increasingly dissipated by a Coulomb-like friction mechanism 
presumably at the gear-teeth interface. Due to this increasing energy loss as the velocity 
nears the resonant areas, fluctuations in velocity increase but the average velocity remains 
roughly the same. Because of this self-reinforcing frictional mechanism, the harmonic- 
drive velocity is prolonged at the frequency at which kinematic error fluctuations excite 
system resonance. This problem is heightened by the non-linear nature of the stiffness 
since increases in torque fluctuation as the system approaches resonance also increase the 
natural frequency of vibration. This effect further prevents the input current from pushing 
through the harmonic-drive resonances since the temporal resonant frequency increases 
alongside the increase in the rotational velocity and resulting spatial excitation frequency, of 
the harmonic drive. In summary, because of the escalating frictional losses and the 
increase in natural frequency which occur as rotational velocities approach system 
resonance, a wide range of input current steps can fall victim to these resonance traps. As 
an explanatory note, the velocity truncation present in the response of the highest current 
step on each velocity plot is due to the saturation limits of the current amplifiers and should 
not be attributed to resonance losses. 

In partnership with the frictional losses which occur at system resonance, 
unpredictable jumps in rotational velocity frequently appear when velocities exit resonant 
areas. As illustrated for joint 1 in figure D.1.2, the step response for a 4.0-amp current 
command holds a relatively constant velocity of about 8000 degrees/second at resonance 
for two seconds then unexpectedly jumps to a final velocity of about 15000 
degrees/second. Additionally, for joint 3, figure D.3.2 illustrates similar behavior for the 
current steps of 0.28, 0.48, and 0.52 amps. This unfortunate dynamic behavior can be 
blamed on three familiar physical mechanisms: (1) the non-linear characteristic of the 
harmonic-drive stiffness profile, (2) the increase in energy dissipation which occurs at 
system resonance, and (3) the large variation in friction which occurs over one output 
revolution of the harmonic drive. First, due to the stiffening-spring characteristic of the 
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harmonic drive, a resonance singularity often appears when the frequency of resonance 
excitation slightly exceeds the natural frequency of vibration. More specifically, as the 
amplitude of torque vibration increases, the equivalent stiffness of the harmonic drive also 
increases as does the natural frequency of system vibration resulting in a frequency 
spectrum similar to the one show in figure 2.45. However, when the harmonic-drive 
velocity breaks through this resonant region, resonance vibration drops immediately due to 
the unstable decline in both torque-vibration amplitude and natural frequency of resonance. 
This non-linear dynamic behavior is commonly referred to as a jump resonance. The 
second accomplice in this jump-resonance behavior is the frictional loss at resonance. Due 
to these losses, as discussed in the previous paragraph, velocity is restricted in resonant 
regions. However, when the harmonic-drive jumps to non-resonating behavior, the 
resonance losses disappear and the leftover energy increases rotational velocity 
dramatically. If the non-linear stiffness and resonance losses were the only factors 
influencing jump-resonance behavior, it is likely that this unpredictable response would 
only appear a very small range of input current commands. However, as observed in the 
experimental results, these jump-resonant responses occur for a surprisingly large range of 
input currents and appear far too often during normal harmonic-drive operation. The origin 
of this unfortunate consequence can be found in the large variation in frictional resistance 
that can exist in harmonic drives over one output rotation. Due to this variation in energy 
dissipation, all operating velocities somewhat near jump-resonance behavior can harness 
enough energy to break through the resonance when frictional variations with output 



Figure 2.45: Non-linear resonance frequency spectrum 
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rotation decrease. Depending on the magnitude of these frictional variations, a wide range 
of operating velocities can be kicked out of resonance due to this decrease in friction. In 
conclusion, due to the interaction of non-linear compliance, resonance losses, and friction 
variation, the stability of harmonic drive operation can be seriously undercut by drastic and 
sudden jumps in operating velocity. 

As noted above to be a contributing agent in jump-resonance behavior, the 
substantial variation in harmonic-drive friction with output rotation can also be observed in 
the velocity time-response plots. As illustrated in all of the response curves of figure 
D.2.2, the variation observed on joint 2 can cause tremendous low-frequency velocity 
variation over one output revolution. To a lesser but noticeable extent, similar variations 
also appear in the response curves for joints 1 and 3 in figures D.1.2 and D.3.2. These 
velocity fluctuations were used to characterize the value of the friction-torque variation on 
each harmonic-drive specimen as summarized in table 2.12. As the low-frequency 
oscillations in velocity indicate, this frictional variation can be expectedly non-sinusoidal. 
The most likely source of this unfortunate behavior is variations in the tightness of fit 
between the flexspline and circular-spline gear teeth, which can be caused by several 
factors such as improper mounting, component misalignment, and manufacturing errors. 


Table 2.12: Variation i 

n output friction torque with output angle 


—sw 

1 ■■■ipii iiiiiiiiiiiiiii 

Torque half-amplitude (N-m) : 

0.6 

1.3 0.03 


2.6.2.1.2 Influence of Observed Behavior on System Properties 

Based on the understanding of typical step-response characteristics gained from the 
input-velocity response, the behavior of the measured harmonic-drive position, current, 
torque, position error, and output velocity can now be better understood. In particular, as 
illustrated in figures D.1.5, D.2.5, D.3.5, D.1.6, D.2.6, and D.3.6, the harmonic-drive 
output velocity response behaves very similarly to the input velocity except that vibrations 
due to kinematic-error fluctuations and system resonance are much larger. Note that the 
joint 3 output velocity is negative since the flexspline, rather than the circular spline, is used 
for output rotation on that joint. Additionally, on the joint 3 output-velocity response, 
distinct periodic velocity fluctuations appear at higher velocities. These fluctuations are 
caused by the output sensor pinion and should not be attributed to harmonic-drive 
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behavior. Also note that the beating in output velocity oscillation amplitude at high 
velocities seen on joints 1 and 2 is due to aliasing in the sampled data rather than an actual 
physical mechanism. 

The current response measured from the motor amplifiers is illustrated in figures 
D.1.3, D.2.3, and D.3.3. Since the motor current is linearly proportional to the motor 
torque, as dictated by the motor torque constant, the behavior of the harmonic drive input 
torque can be understood. From the response curves for each joint, it can be seen that at 
low current commands, the amplifiers can sustain a fairly constant input torque in the face 
of harmonic-drive torque fluctuations. However, at higher current commands, the current 
amps have difficulty keeping pace with the higher frequency torque fluctuations and 
noticeable variations in current amplitude occur. The highest current-command response in 
each figures illustrate the vulnerability of the current signal when the amplifier saturates 
from the high motor velocity. Due to imperfections in the A to D board, significant 
electrical noise appears in the current signals on joints 2 and 3. Fortunately, this noise had 
no effect on the dynamic response of the harmonic-drive systems. 

The output torque response measured from the torque sensor on each joint is 
illustrated in figures D.1.7, D.2.7, and D.3.7. As a reminder, because of the location of 
the torque sensor on each joint, this measured response is a direct reading of the torque 
experienced at the flexspline of the harmonic drives. This output torque has two 
contributions: (1) the dynamic torques generated by the acceleration of the output inertia, 
and (2) the damping torque produced by the output bearings. Since the damping torque in 
the output bearings is very small compared to the inertial torques generated by the position- 
error fluctuations, the measured torque response in all cases oscillated around zero at all 
velocities. As shown on all three joints, torque fluctuations can become extremely large at 
velocities which excite system resonance. For example, on the largest harmonic-drive 
system, joint 1, torque readings reached a powerful 80 N-m during system resonance. As 
expected, the frequency of torque fluctuation is identical to the frequency of velocity 
fluctuation which is equal to the rotational speed of the input times the position-error 
excitation frequency at the appropriate resonance. 

The input and output position response measured on the three harmonic-drive 
specimens is illustrated in figures D.1.1, D.2.1, and D.3.1, and figures D.1.4, D.2.4, and 
D.3.4. From these plots, it can be seen that substantial fluctuation noticed on the velocity 
plots is virtually unnoticeable on position measurements. However, due to increased 
frictional losses at resonances, the increases in final position for constant increases in the 
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current-step input is by no means constant. For all tests, the input and output position 
responses display identical behavior except that they differ in magnitude by the given 
harmonic-drive reduction ratio. 

As shown in figures D.1.8, D.2.8, and D.3.8, the position-error response of a 
single trial is shown for the three harmonic drives. The given plots were selected to 
illustrate typical position-error behavior near a system resonance and were derived by 
subtracting the output position from the input position divided by the gear ratio for the 
given step-response trial. The important feature to note in these plots is the amplification of 
position error at system resonance due to torsion across the harmonic-drive transmission. 
As the plots illustrate, this effect can increase the dynamic position inaccuracy measured 
across the drives by an order of magnitude over static tests. This large effect of dynamic 
torque fluctuations on the measured position-error illustrates the importance of static tests to 
ensure accurate measurements of harmonic-drive kinematic error. As an explanatory note, 
the low-frequency oscillations seen in the position-error response are due to manufacturing 
errors in the output sensor gear and pinion and should not be attributed to the harmonic 
drive. 


In general, changes in the dynamic response for identical step-response trials was 
minimal as long as harmonic-drive parameters remained unchanged. For example, since, 
as discussed earlier, factors such as harmonic-drive preload, assembly, lubrication, input 
and output orientation, and operating temperature, can greatly influence system parameters 
such as stiffness and frictional losses, changes in any of these conditions between trials can 
greatly influence dynamic response. Care was taken to make sure that these parameters 
remained consistent across all step-response trials on each joint. 


2.6.2.2 Modeling Considerations 

Based on the observations made about typical harmonic-drive step-response 
behavior, theories can be developed about how to model the observed dynamic effects. 
First, in order to accurately predict resonance behavior, the input and output inertias, 
kinematic error-signature, and the harmonic-drive stiffness profile must be carefully 
measured and introduced into a dynamic model. Second, as described above, to capture 
the observed jump-resonance behavior, three physical mechanisms must be modeled: (1) 
the non-linear stiffness profile, (2) the Coulomb-like frictional losses at resonance, and (3) 
the variation in friction with output rotation. The first and last of these three items can be 
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measured relatively easily and modeled by mathematical functions which approximate the 
observed behavior. The second item, however, relies on the interaction between torque 
fluctuations in the harmonic-drive and rubbing at the gear-tooth interface. Therefore, any 
model that is supposed to accurately describe frictional losses at resonance must include a 
representation which predicts the dynamic forces acting on the gear-tooth geometry. 
Additionally, as illustrated in the motor-current response curves, a dynamic model of the 
current amplifiers might be necessary in some domains to capture the saturation behavior 
and dynamic fluctuations. 


2.6.3 Dynamic Response Conclusions 

The unusual and surprising dynamic behavior observed in harmonic drives is the 
result of the interaction of several transmission mechanisms. First, due to the kinematic 
error in the harmonic drive, resulting dynamic torque fluctuations amplified by the 
transmission ratio can be substantial during operation. When these large torque 
fluctuations are combined with the low stiffness of the transmission, significant resonance 
vibration can occur in typical operating ranges. Since coulomb friction at the gear-teeth 
interface is amplified by dynamic torque fluctuations, extra energy is dissipated at 
resonance and operating velocities are suppressed for increasingly larger motor currents. 
However, when the input motors supply enough energy to break through these resonant 
regions, a sudden reduction in resonance losses results in a sharp increase in velocity and a 
significant decrease in torque fluctuations. Since the stiffness profile in harmonic-drives is 
non-linear, the natural frequency of the system increases with increased torque levels. This 
non-linear effect has an adverse affect on dynamic behavior around system resonance since 
the increase in natural frequency prolongs operation at resonance and accentuates the 
reduction in vibration when these resonant regions are surpassed. Lastly, since harmonic 
drives can have a large variation in friction torque over one output revolution, the resulting 
substantial variation in operating velocity can greatly influence the location of the jump 
resonance. In order to accurately model the dynamic behavior of the harmonic drive, 
appropriate representations and interactions of these individual effects must be understood. 
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The goal of modeling any physical mechanism is to discover the simplest 
representation which can replicate system performance to the desired level of accuracy. 
When modeling a mechanical assembly with a harmonic-drive transmission, several levels 
of detail can be incorporated into the harmonic-drive model to improve the integrity of the 
system representation. In particular, as seen in the experimental results, harmonic-drive 
properties which can significantly influence system performance include 



static, dynamic, and periodic frictional losses, 


1 

1 


non-linear transmission compliance, and 
kinematic transmission error. 


The goal of this section is to characterize the modeling accuracy that can be achieved by 
incorporating different combinations of these transmission properties into a harmonic-drive 
representation. From the complex behavior observed in the measured dynamic-response 
data, finding a model which accurately characterizes harmonic-drive behavior proves to be 
a formidable task. However, if such a representation can be created, design and control of 
systems containing harmonic drives can be substantially improved. Additionally, the 
knowledge embedded in such a dynamic model can become a powerful medium for 
understanding and improving the design and operation of harmonic-drive transmissions. I 
hope that the results presented in this section can equip and enlighten the expanding 
community of harmonic-drive users and designers. 

My modeling efforts have focused on developing a series of harmonic-drive models 
at different levels of complexity and evaluating their performance. First, I analyzed a 
transmission model that ignored all non-ideal transmission properties and discovered that it 
provided an extremely poor description of actual dynamic performance. Second, model 
accuracy improved when frictional losses were included into the representation, but the 
torque and velocity fluctuations observed in the experimental response did not appear in the 
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theoretical predictions. To address this issue, compliance and then kinematic error were 
included into the transmission model. This improved representation produced simulated 
results which more closely matched the experimental observations but failed to capture the 
significant energy losses which occur at system resonance. To eliminate this discrepancy, 
a final model which incorporated gear-tooth geometry and coulomb friction was developed 
which improved predicted performance around resonance. Unfortunately, due to the 
sensitivity of this model to harmonic-drive properties and the inherent difficulty in 
accurately measuring these properties experimentally, further research must be done to 
refine this representation into a viable prediction tool. Nevertheless, the novel framework 
developed in this model provides a simple and powerful representation for hatching and 
nurturing an intuitive understanding of harmonic-drive operation. 

A complete description of my theoretical analysis of the harmonic-drive 
transmissions is presented in this section. To establish the foundation of my conceptual 
investigation, I will first explain the notation conventions used throughout this section and 
describe the numerical simulation used the dynamic equations. Next, the modeling 
representations used for the three harmonic-drive testing stations will be discussed apart 
from the individual transmission models. In light of this framework, I will then develop a 
series of dynamic models of increasing complexity to illustrate the range of performance 
delivered by different harmonic-drive representations. For each of these models, equations 
of motion will be derived, operating parameters will be calculated from experimental 
measurements, and the accuracy of simulation results will be characterized through 
comparison with actual harmonic-drive behavior. On the basis of these theoretical 
observations, conclusions and recommendations will be made for each harmonic-drive 
model about its usefulness and range of applicability. 


96 


Chapter 3: Modeling Techniques and Noiation Conventions 


3.1 Modeling Techniques and Notation 
Conventions 

Modeling techniques and conventions will remain consistent for all of the dynamic 
models presented in this section. Specifically, Newton’s method will be employed 
exclusively to determine dynamic force interactions and derive equations of motion. The 
position, velocity and torque variables used in these equations will be illustrated in a figure 
for each model and governed by the following conventions: 

0-variables represent rotational position, 
co-variables indicate angular velocity. 

T-variables denote torque. 

All variables are referenced from ground. 

All positions and velocities are defined to be positive in the clockwise 
direction when viewed from the input-side of the harmonic drive. 

All torques are defined to be positive in the direction in which they will act 

when a positive load is applied to the input inertia and the output inertia 
remains stationary. 

The positions, velocities, and torques seen by the three ports of the harmonic 

drive will have the subscripts wg for the wave-generator, fs for the flexspline, 
and cs for the circular spline. 

Additional variables such as model parameters and constants that are used for the testing- 
station models are defined in table 3.1. Similarly, table 3.2 lists descriptions of the 
variables which will appear in the five different harmonic-drive models to be developed 
later. Sporting this notation lifejacket, I will now plunge into the depths of equation 
derivation and analysis for five different harmonic-drive models. 

In order to accommodate different harmonic-drive models, the equations of motion 
for the testing-stations will treat the transmission as a black box with three interacting ports: 
the wave-generator, the flexspline, and the circular-spline. Velocities and torques will be 
transmitted through these three ports but the testing station models will be blind to the 
actual dynamics inside the transmission. Following the complete derivation of the testing- 


0 

0 

0 

0 

0 


2 


97 



Chapter 3: Modeling Techniques and Notation Conventions 


Table 3.1: Explanation of variable names for testing-station models 


i m current sent to the DC motors by the amplifiers 

J in input inertia (motor armature and wave-generator) 

J ou t output inertia (load mounted at harmonic-drive output) 
bi n linear damping coefficient for input friction 

bout linear damping coefficient for output friction 

K t motor torque constant 

K b motor back-EMF constant 


Table 

3.2: Explanation of variable names for harmonic-drive models 

N 

harmonic-drive catalog gear-ratio 

Tb_constant 

velocity-independent friction torque in the harmonic-drive 

Tb_dynamic 

velocity-dependent friction torque in the harmonic-drive 

T b_cyclic 

friction torque in the harmonic-drive which varies with output position 

bconstant 

value of the velocity-independent friction torque 

bl 

linear damping coefficient 

b2 

cubic damping coefficient 

A b 

half-amplitude of the cyclic-friction torque function 

<l>b 

phase of the cyclic-friction torque function 

ki 

linear stiffness coefficient 

k 2 

cubic stiffness coefficient 

Al 

half-amplitude of the first kinematic-error component at output 

A2 

half-amplitude of the second kinematic-error component at output 

A 3 

half-amplitude of the third kinematic-error component at output 

01 

phase of the first kinematic-error sinusoid 

02 

phase of the second kinematic-error sinusoid 

03 

phase of the third kinematic-error sinusoid 


Coulomb friction coefficient at the harmonic-drive gear-tooth interface 


station models, five different harmonic-drive models will be developed that can be inserted 
easily into the overall system models. Since the inertias of the wave-generators are 
included with the motor armature inertias, no inertia elements are needed in any of these 
harmonic-drive representations. Consequently, the equations of motion for each will 
consist solely of algebraic rather than differential equations. These algebraic equations will 
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outline the relationships between the torques and rotations on the three harmonic-drive 
ports. In particular, since the inertias in the testing stations dictate the position and velocity 
seen at the input and output of the harmonic-drive, equations will be developed to relate 
these rotations to the torques experienced by the transmission. Additionally, because these 
equations completely characterize the three-port behavior of the harmonic-drive, any two- 
port configuration of the transmission can be represented easily by setting the rotation of 
the non-rotating harmonic-drive port to zero. 

Due to the complexity of the differential equations generated when the harmonic- 
drive models are inserted into the overall system models, analytical solutions describing 
dynamic behavior were rarely available. However, using a Runga-Kutta numerical 
integration scheme with adaptive step-size control, the dynamic behavior specified by a set 
of differential equations could be simulated over a given time period for a specific set of 
input conditions. The code which was developed to generate dynamic time-response data 
for different harmonic-drive models is listed in Appendix B. In addition to calculating the 
forces, positions, and velocities predicted by the differential equations, this simulation code 
determines the energy distribution in the system and can be used generate static stiffness 
data as well. 

To facilitate the development of new harmonic-drive representations in the computer 
simulation, the modularity of the different harmonic-drive models was preserved by 
collecting the equations for each in a single function. These harmonic-drive functions can 
be used to calculate the torque on the three harmonic-drive ports when the corresponding 
positions and velocities are specified. If the proper input parameters and variable 
definitions are provided, these harmonic-drive functions can be extracted easily from the 
code in Appendix B and used in other dynamic simulations. 

When executed on a Sun SPARCstation ELC, typical simulation times for a four- 
second time-response using the more complex harmonic-drive models were much less than 
one minute. Although the exact time for simulation is influenced not only by the 
complexity of the dynamic equations but by the desired accuracy and nature of the system 
parameters as well, in general, the small computational requirements of the harmonic-drive 
models presented below suit them well for rapid numerical simulation. Additionally, it is 
likely that simulation performance can be further improved through refinement of the code 
which calculates values using the dynamic equations. 
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3.2 Modeling the Experimental Apparatus 

To distinguish the dynamic effects of the harmonic-drive transmission from 
dynamic behavior observed in the experimental apparatus, dynamic models of the testing 
stations are required. Since, of the three testing assemblies, two different design 
configurations exist, two separate dynamic models must be constructed. Once this is done, 
different harmonic-drive representations can be inserted into these models, and dynamic 
behavior of the aggregate system can be simulated for a range of transmission models. 
This section is devoted to describing the two models developed to capture the behavior of 
the experimental equipment and deriving the corresponding equations of motion. 


3.2.1 Description of Models 

As described in section 2.1, two different designs were incorporated into the 
harmonic-drive testing facilities in order to observe harmonic-drive behavior in different 
operating configurations. On joints 1 and 2, the wave-generator is driven by the input 
motor, the circular spline carries the output rotations, and the flexspline is mounted to 
ground. On joint 3, the wave-generator is still driven by the input motor, but the output is 
conveyed through the harmonic-drive flexspline while the circular spline is attached to 
ground. These two different joint designs are illustrated in figures 3.1 and 3.2. 

The relative location of important model parameters can be understood from these 
figures. First, for both joint configurations, the DC motor applies a torque to the input 
inertia which is composed of the motor armature and the harmonic-drive wave-generator. 
Additionally, each joint sports an output inertia consisting of a large mass mounted to the 
end of the output link on the output port of the transmission. Due to motor damping, a 
friction torque can be seen between the motor shaft and the motor housing, and similarly, 
the output bearings introduce a finite friction torque between output rotation and ground. 
Since, the motor housing on joints 1 and 2 rotates with the circular spline, the absolute 
rotation seen by the input inertia is equal to the sum of the input and output rotation. The 
input inertia on joint 3, on the other hand, sees the direct rotation of the motor shaft relative 
to ground. 
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Based on this understanding of the motion constraints and operating parameters, 
two simple lumped parameter models distilled from the two joint designs are shown in 
figures 3.3 and 3.4. As illustrated in these figures, the harmonic-drive model is 
represented as a black box with three ports that transmit the torque and rotation seen by the 
wave-generator, flexspline, and circular spline. Different models that capture harmonic- 
drive behavior in varied levels of detail will be developed later and substituted into this 
black box. 
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3.2.2 Derivation of Equations of Motion 


Using the lumped-parameter models illustrated and described above, the generation 
of dynamic equations is straightforward. The variable names shown for both models in 
figures 3.5 and 3.6 will be used throughout this derivation and governed by the notation 
conventions itemized in section 3.1. For further clarity, important model parameters were 
summarized previously in table 3.1. Using the given notation, naming conventions, and 
parameter definitions, the equations of motion for joints 1 and 2 will now be derived 
followed by a similar derivation for joint 3. 
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By applying Newton’s Law to the input and output inertias in the model shown in 
figure 3.3, the resulting equations of motion for joints 1 and 2 are 

Jin®in = T m - - T W g , (3.1) 

and 


Jout®out — T cs - T b_ out + T b in • (3.2) 

In order to solve these equations, the torques on the right-hand side of each equation must 
be expressed in terms of the state variables, 0 in and 0 out . This can be done readily for the 
damping torques by assuming that they are linear functions of velocity: 
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Tb_in — bi n (Bin " Bout) » &nd 

(3.3) 

Tb_out =: bout Bout • 

(3.4) 


Since the motor torque is linearly proportional to the current produced by the motor 
amplifiers, T m can also be determined: 

T m = K l 1 m, (3.5) 

where the amplifier current, i m , can be found from the equations presented in Appendix C 
which describe the saturation limits of the amplifiers. The remaining unknown torque 
components, T wg and T cs , can be found by specifying the kinematic constraints imposed by 
the system on the three harmonic-drive pons and calculating the resulting torques using one 
of the harmonic-drive models presented later in this section. Since each port of the 
harmonic drive is directly coupled to either ground or one of the inenias, these kinematic 
constraints are simply 


0 W g — Bin , 

(3.6) 

0f s = 0 , and 

(3.7) 

0CS = Bout • 

(3.8) 


Therefore, given these equations, and any dynamically complete harmonic-drive 
representation, the equations of motion are fully defined. 

The equations of motion for the joint 3 testing station are very similar to those listed 
above for joints 1 and 2. Using figure 3.4, they can be determined readily to be 

Jin^in = T m - Tbj n - T wg , and (3.9) 

JoutQout = ' Tf s - Tb_ ou t . (3.10) 

The output damping torque for this model can still be calculated using equation (3.4) 
above, but the input damping torque must be determined from 

T'b_in = bj n Bin • (3.1 1) 
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Additionally, the motor torque expression remains the same as equation (3.5) above, but 
the kinematic constraints on the harmonic-drive ports are now 


0\vg — ®in > 

(3.12) 

0fs = Oout i and 

(3.13) 

e cs = o. 

(3.14) 


Therefore, by substituting these constraints into a harmonic-drive model, calculating the 
resulting torques on the harmonic-drive ports, and substituting them into the equations of 
motion listed above, the dynamic behavior of the joint 3 testing station is fully defined. 


3.2.3 Model Parameter Calculation 

Having derived the governing equations for the two different harmonic-drive testing 
stations, dynamic response can be simulated once values have been determined for the 
given model parameters. First, the input inertia can be found by summing the catalog 
inertia values for the wave-generator and the motor armature, and the output inertia can be 
calculated from the weight and location of the output mass on each joint. Since all other 
output inertia components on each joint were negligible compared to the inertia of the 
output mass, they were ignored in the output inertia calculations. All important inertia 
values are summarized in table 2.10 and reproduced in table 3.3 below. Second, a value 
for the input damping coefficient can be obtained from the DC motor catalog, and the 
damping in the output bearings can be found from roller-bearing catalog specifications. 


Table 3.3: Paramet 

er values for the 

testing-station n 

nodels 

Jin 

2.2e-4 

joint 2 

7.6e-5 

joint 3 

5.2e-6 

Jntt (kB-rn 2 ) 

0.268 

0.056 

0.0088 

fci rt (N-m/{lnput deg/sec)) 

2.3e-7 

2.3e-7 

1 .Oe-7 

&<wf {N-ro/(output dog/sec)) 

<7.0e-4 

< 7.0e-4 

< 4.0e-4 

K* (N-m/amp) 

0.06 

0.03 

0.04 

Kfe <voKs/{tnput deg/sec)) 

1.05e-3 

5.07e-4 

7.08e-4 
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These values are also listed in table 3.3. Lastly, as mentioned above, the torque produced 
by the DC motor is linearly proportional to the current it receives from the amplifiers. This 
current is specified by the desired input command and remains true to this value unless high 
motor velocities cause the current amplifiers to saturate. Using the catalog motor 
specifications listed in table 3.3 and the measured motor velocity, appendix C presents a 
simple algorithm which can be used to capture saturation behavior. Having developed 
representations which describe testing-station operation and determined values for 
corresponding model parameters, attention can now be focused directly on the harmonic- 
drive transmission. 
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3.3 Harmonic-Drive Model 1: 
Ideal Transmission 


At a bare minimum, a harmonic-drive model should be able to describe the overall 
gear-reduction properties of the transmission. Since the harmonic-drive has three 
independent rotational components, equations which describe three-port transmission 
behavior must be used. Based on the information provided in manufacturer’s catalogs, a 
simple three-port harmonic-drive transmission model, as illustrated in figure 3.7, can be 
characterized by the following equations: 

0 wg = (N+ l)0 C s-N0 fs , (3.15) 


co wg = (N + 1) co cs - N (i)f s , and 


(3.16) 


T = 

A \i /a — 


wg (N+l) 


T cs = tr T fs • 


(3.17) 


Note that a positive value of Tf s will cause a negative rotation of the flexspline since it is 
defined to be positive in the direction opposite to T wg and T cs . As previously discussed, N 
is the catalog gear ratio and the remaining variables are illustrated in figure 3.8. Using 
these equations, the positions, velocities, and torques for any harmonic-drive configuration 
are completely defined. 


When this representation is substituted into the different joint-configuration models 
presented above, the input inertia becomes directly coupled to the output inertia creating a 
first-order system. Consequently, any step input in motor current should produce an 
exponential increase in velocity up to the point where the motor torque matches the dynamic 
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damping torque of the motor and output bearings. This behavior is illustrated in the 
simulated input velocity responses for the three harmonic-drive testing stations shown in 
figures 3.9, 3.10 and 3.11. 


By comparing these predicted input-velocity responses to the experimental input- 
velocity responses shown in Appendix-D-figures D.1.2, D.2.2, and D.3.2, it can be seen 
that the steady-state velocities achieved using the ideal transmission model are several 
orders-of-magnitude higher than the actual system. Additionally, the velocity fluctuations 
observed in the experimental results do not appear in the model predictions. These two 
unsurprising discrepancies can be accounted for by the lack of any frictional losses or 



110 






Chapter 3: Harmonic-Drive Model 1 






Chapter 3: Harmonic-Drive Model 1 


kinematic error in the transmission model. Given this justification for poor model 
performance, it can be concluded that the only circumstances where this representation may 
prove useful is when harmonic-drive frictional losses are small compared to other losses in 
the mechanical system and when velocity and torque fluctuations are unimportant. For 
example, under static loading, when dynamic friction is zero, this model can be used to 
accurately predict torques when applied loads are significantly greater than the static friction 
in the transmission. Additionally, in static positioning applications where the desired 
accuracy is lower than the kinematic error in the transmission, this model may be useful to 
calculate the rotational displacement of the harmonic-drive ports. 
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3.4 Harmonic-Drive Model 2: 

Transmission with Frictional Losses 

As concluded from the poor results of the ideal model, energy losses in the 
transmission must be acknowledged in order resurrect the simulated velocity response. By 
matching the frictional torque in the model to the motor torque resulting from the input 
current, accurate recreations of experimental steady-state velocity values can be achieved. 
Furthermore, by including a periodic friction function which varies with the position of the 
transmission output, the undulations in steady-state velocity observed in the experimental 
results can also be replicated. As discussed below, by incorporating these simple 
modifications, the response of the ideal transmission model can be improved substantially. 

3.4.1 Description of Model 2 

The bulk of energy dissipation in harmonic-dnves can be blamed on two sources: 
(1) losses in the wave-generator bearing and (2) friction due to gear-tooth meshing. Of 
these two candidates, it is probable that the majority of the frictional dissipation results 
from gear-tooth-rubbing action between the flexspline and circular spline. Based on this 
assumption, a new harmonic-dnve model is shown in figure 3.12 in which a dissipation 
element is included between the flexspline and circular spline. Since the location of these 
frictional losses in the model can greatly influence simulated behavior, careful investigation 
was performed to justify this representation. 
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In addition to conceptually capturing the rubbing losses between the flexspline and 
circular spine gear teeth, the location for the frictional losses shown in figure 3.12 was 
selected for several additional reasons. First, as shown in figures D.1.7, D.2.7, and 
D.3.7, since torque-sensor readings in experimental dynamic response measurements 
display very large fluctuations but virtually no frictional bias torque, it can be concluded 
that the torque sensors are blind to any frictional torque present in the harmonic drive. On 
joint 1, for example, if the friction losses in the model were placed between the circular 
spline and ground, the torque sensor mounted to the flexspline would not only see any 
torque fluctuations on the output inertia but would also directly measure the frictional 
torque from the transmission losses. By moving these frictional losses to act between the 
flexspline and circular-spline, the torque sensors for both joint configurations ignore these 
losses and return results that more closely resemble the experimental observations. As 
another justification for the selected location of transmission losses, any static friction 
included in this dissipation element will be manifested as hysteresis loss in a stiffness 
profile collected from either output port of the transmission. On the other hand, if the 
friction mechanism occurred between the wave-generator and either output port, hysteresis 
loss would appear only in the stiffness profile collected from the output port connected to 
that friction element. From this distinction, it seems highly probable that frictional torques 
in the harmonic-drive occur primarily between the flexspline and circular spline gear-teeth 
and can account for the hysteresis loss observed in experimental stiffness profiles. To 
substantiate this claim empirically, results were gathered from models employing different 
friction mechanisms and the integrity of the representation derived in this section was 
verified. 


3.4.2 Derivation of Dynamic Equations 

The addition of frictional losses to the ideal transmission model complicates the 
dynamic equations only slightly. Since, using this new representation, the torque seen by 
the flexspline and circular spline is the combined torque of the ideal transmission element 
and the frictional torque, the following equations result: 

Tf s = T n _f s - T b , and (3.18) 

Tcs = T n cs - Tb . ( 3 . 19 ) 
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T b =f(co) 


Figure 3.13: Notation tor the harmonic-drive model with friction 



In these equations, Tb is the aggregate friction-torque function and, as implied by the sign 
conventions stated previously, T n _f s is defined to be positive in the opposite direction of 
T n _cs- The variables used in these equations and the rest of this derivation are illustrated in 
figure 3.13. Using the ideal-transmission torque relationship in equation (3.17), the torque 
on the wave-generator can be directly related to T n _f s and T ncs : 



1 

(N+l) 


n cs 



(3.20) 


Given these three equations, the torque constraints on the model are complete, 
however, the friction-torque function, Tb, must be defined in order to finish the 
representation. As observed in experimental results, harmonic-drive friction typically has 
three components: (1) a velocity-independent friction term, (2) a non-linear velocity - 
dependent damping function, and (3) a periodic friction torque which varies with output 
rotation. Experimental observations also suggest that frictional losses can be greatly 
influenced by other parameters such as preload, wave-generator orientation, and operating 
temperature. For the sake of model simplicity, these frictional influences will be ignored. 
Additionally, experimental results show that harmonic-drives exhibit a substantial static- 
friction torque which must be overcome before rotation can occur. Since this stiction 
torque is also influenced by several parameters such as preload and input angle, an accurate 
model of this property relies heavily on careful experimental characterization. Once this 
experimental understanding is complete, a stiction model can be implemented by placing 
appropriate constraints on the minimum motor torque required to effect rotation. Although 
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this model was developed in the computer simulation, it will be ignored here since it only 
influences dynamic behavior when velocities are close to zero and consequently has little 
effect on the experimental and theoretical dynamic response. Given these assumptions, the 
friction function which provided optimal results with minimal complexity is defined as 

Tb — Tb_ cons ta n t b Tb_dynamic b Tb_cyclic > (3 2 \) 


where 


Tb_constant ~ h CO nstam ) 

Tb_dynamic = bj (Cl) n _ cs - ^n_fs) b t>2 (C0 n _ _cs ~ CDn fs) 3 , and 
Tb_cyclic — Ab sin ((0n_cs “ 0n_fs) b (j>b) . 


(3.22) 

(3.23) 

(3.24) 


As equation (3.21) states, this model captures the three facets of friction behavior noticed in 
the experimental results. First, the velocity-independent friction component, as described 
in equation (3.22), is represented by a single value, b constam . Second, since experimental 
data indicated that velocity-dependent friction could be accurately captured by a cubic 
relationship, the dynamic-friction torque defined in equation (3.23) is a cubic function of 
the relative velocity between the flexspline and circular spline. In this equation bi and b 2 
represent the linear and cubic damping coefficients, respectively. Since all observed 
experimental behavior indicates that dynamic damping is highly non-linear, a linear 
damping representation is likely to cause inaccurate model performance. Third, periodic 
variations in the frictional torque are captured in the sinusoidal function in equation (3.24). 
As this relationship indicates, frictional torque fluctuations of half-amplitude Ab complete 
one cycle every time the circular-spline makes one complete rotation relative to the 
flexspline. To match this model to experimental observations, a phase-shift of <j)b is also 
included in this equation. Using the experimental results from sections 2.4, 2.5, and 2.6 
of this document, values for the friction parameters b CO nstant> b], b 2 , A b , and 0 b can be 
determined readily. 


With the given definition of the harmonic-drive friction, the transmission model is 
virtually complete. The only remaining task is to realize that the variables 0 n _f s , 0 n _ C s> 
co n _f s and co n _ cs are identical to the variables 0 fs , 0 CS , co fs and co cs respectively, due to the 
kinematic constraints in the model. Given this identity, the friction-torque function can be 
calculated based on the rotations specified on the harmonic-drive ports and the remaining 
torques in the transmission can be found. 
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3.4.3 Evaluation of Simulated Results 

Using values for the harmonic-drive friction parameters taken directly from the 
experimental observations, several numerical simulations were performed for each 
harmonic-drive over a range of current commands. As a sole indicator of model 
performance, the resulting input-velocity response for the three different harmonic drives is 
illustrated in figures 3.14, 3.15, and 3.16. These three plots can be compared to the 
experimental input-velocity response shown in figures D.1.2, D.2.2, and D.3.2 in 
Appendix D to evaluate model integrity. 

By comparing the theoretical and experimental time-response, it is clear that the 
harmonic-drive model which includes frictional losses boasts performance greatly 
improved over the ideal transmission representation. In particular, for high velocities, 
input velocity predictions are almost identical to experimental observations on each joint. 
For all these curves, response is very similar to first order-system response where the time- 
constant is dictated by the total system inertia and the steady-state velocity is reached when 
the motor torque matches the total friction torque. The only discrepancies between these 
high-velocity response curves is in the shape of the velocity fluctuations due to friction 
variations with output rotation. By replacing the simple sinusoidal friction in the model 
with a more complex function which closely resembles the actual torque variations in the 
transmission, it is likely that theoretical predictions could replicate the actual response 
almost exactly. Low-velocity response predictions also show promise since theoretical 
steady-state velocities and curve shapes show good agreement with experimental results. 

Despite these modeling successes, simulation performance plummets around 
velocities which excite system resonance. Specifically, the substantial velocity fluctuations 
which pervade the experimental velocity response are absent in the theoretical predictions. 
These velocity fluctuations can become extreme around resonance and utterly dominate 
system behavior. Since the model fails to capture these disturbances, steady-state velocities 
and curve-shapes are markedly different from actual results. Additionally, as illustrated in 
figures D.1.5, D.2.5, and D.3.5, the more prominent fluctuations in output velocity are 
substantial even at high velocities. Since the output-velocity response predicted by the 
model, although not shown here, will capture none of this agitated behavior, simulation 
accuracy struggles through all ranges of operation. As a final criticism of the virtue of this 
harmonic-drive model, the periodic velocity fluctuations due to cyclic friction, as shown in 
figure 3.15 in particular, are smaller in amplitude at high velocities than the experimental 
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results. However, at low velocities, these periodic variations are in closer agreement with 
experimental response. This unfortunate behavior indicates that the amplitude of the cyclic 
friction function is not constant as anticipated, but may indeed vary with such variables as 
rotational velocity. 


3.4.4 Model 2 Conclusions 

Based on the performance of the harmonic-drive model with frictional losses, if 
small velocity fluctuations can be tolerated, it is probable that this model may prove useful 
for predicting the dynamic behavior of harmonic-drive systems designed to operate 
predominantly at velocities far removed from regions of resonant behavior. However, in 
order to ensure the accuracy of this model, a comprehensive experimental friction 
characterization is first required. In particular, transmission friction must be measured for a 
range of velocities and output orientations under the identical environmental and operating 
conditions to be experienced during normal operation. Since the velocity-independent 
friction torque looms relatively large at small rotational velocities, accurate characterization 
of this value ensures good model performance at low velocities. Conversely, dynamic 
friction forces, which become significant at high velocities, govern model performance for 
large current inputs. Lastly, careful measurement of experimental cyclic friction is 
rewarded by accurate reproduction of the resulting velocity fluctuations. By following 
these guidelines, this simple transmission model can be used to make accurate predictions 
for a limited range of harmonic-drive systems. 

In conclusion, the results of this section illustrate that the harmonic-drive model 
with friction performs significantly better than the ideal transmission model. Due to the 
incorporation of observed energy dissipation into the model, average steady-state velocities 
show fairly close agreement for operating regions where resonant vibration is minimal. 
Unfortunately, since the model ignores kinematic position-error variations, theoretical 
predictions are devoid of the velocity fluctuations which appear throughout the 
experimental response. Additionally, because these velocity fluctuations become intense 
during system resonance, theoretical results become very inaccurate over a large portion of 
the observed operating range. Based on these shortcomings, it is clear that an accurate 
model must capture velocity turbulence and resonance behavior. Since, as discussed 
previously, this behavior is spawned from the interaction between kinematic error and 
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transmission compliance, an improved model should acknowledge these two harmonic- 
drive properties. 
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3.5 Harmonic-Drive Model 3: 

Transmission with Friction and 
Compliance 

As published research suggests, once frictional behavior is captured, the first 
strategy adopted by designers to improve performance of harmonic-drive models is to 
include a linear or non-linear compliance in the transmission. Since this compliance 
transforms the first-order model into a second-order system, a vibration mode appears and 
resonance phenomenon can be anticipated. This section outlines this new model and 
demonstrates that, unless large torque loads are present in the representation as well, the 
compliant transmission model does little to improve simulation performance over the model 
with friction alone. 

3.5.1 Description of Model 3 

A lumped-parameter schematic illustrating the harmonic-drive model with 
compliance and friction is illustrated in figure 3.17. Since, as noted in previous research, 
harmonic-drive compliance occurs in the wave-generator as well as the gear-teeth, I chose 
to locate the ideal spring on the input side of the gear-reduction in this model. However, 
by adjusting the coefficients of the stiffness function appropriately, this compliance can be 
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reflected to either of the other two ports of the transmission with no effect on the dynamic 
behavior of the system. Due to the location of the friction element with respect to the 
compliance in this model, the stiffness profile measured from either output port of the 
transmission will display the characteristic hysteresis loss due to the velocity-independent 
friction component in the friction-torque function. 


3.5.2 Derivation of Dynamic Equations 

Given this simple lumped-parameter model and the notation definitions illustrated in 
figure 3.18, the derivation of the equations of motion is straightforward. First, using 
force-balance equations in (3.18) and (3.19) as well as the friction function in (3.21), the 
torque on the output ports of the harmonic-drive, Tf s and T cs , can be related to the variables 

T n _f s and T n _ cs . Second, using the ideal gear-reduction relationship provided in equation 

(3.17), these torques can then be equated to the torque across the spring and the torque at 
the wave-generator: 


T k = T wg 


= T 


n_wg 


=-1_t 

(N + 1) n - cs 



(3.25) 



T b =f(a>) 


Figure 3.18. Notation for the model with friction and compliance 
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Given this identity, the torques on all three harmonic-drive ports can be calculated once the 
torque across the spring is known. 


To calculate the transmission stiffness, several different functions can be used. 
Specifically, depending on the exact shape of the stiffness profile for a particular harmonic- 
drive, a linear stiffness relationship may suffice. However, as observed in the 
experimental stiffness observations, most harmonic-drives, especially under heavy loading, 
display distinctly non-linear compliance. Also as seen in the experimental results, this non¬ 
linear behavior can be captured typically by a cubic relationship: 

Tk = ki (0 W g - 0n_wg) + k2 (0wg " 0n_wg) > (3.26) 

where ki and k 2 are the linear and cubic stiffness constants, respectively. To further refine 
this stiffness characterization, a more complex equation could be used which includes the 
influence of wave-generator orientation and gear-tooth preload on the stiffness profile. 
Since, as experimental observations indicate, this influence can cause fluctuations in 
stiffness profiles of over thirty percent, a truly accurate stiffness function cannot ignore 
these effects. However, since the necessity for this level of detail is unclear, these 
influences will be neglected and the simple cubic relationship in (3.26) will be featured 
exclusively in the dynamic model. 


The derivation of equations for this model can be completed by recognizing that 
kinematic constraints dictate 0 n fs is equal to 0f s and 0 n cs is equal to 0 CS . From these 
identities, 0 n _wg can now be found from the kinematic requirement of the ideal 
transmission model in equation (3.15): 

0n_wg = (N + 1) 0 CS - N 0f s . ^ 

Given values for 0 wg , 0f s , and 0 CS , these equations can be used to calculate the resulting 
torques on the three ports of the harmonic-drive transmission. 


3.5.3 Model 3 Results and Conclusions 

To evaluate the validity of the stiffness representation used in this model, two 
hallmarks can be employed. First, the stiffness function should replicate the experimental 
dynamic behavior observed at system resonance. In particular, resonant frequencies and 
vibration amplitudes in the model should coincide with observed experimental behavior. 
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Second, when static stiffness data is generated from the transmission model, it should 
closely resemble the experimental stiffness profiles. Specifically, if the interaction between 
the model friction and compliance is accurate, the hysteresis loss and curve shape in the 
actual stiffness data should be captured. In order to observe the influence of resonant 
behavior, both of these criteria will be discussed in light of the harmonic-drive model 
which includes kinematic transmission error to be discussed later. 

Using experimental stiffness results to determine the stiffness coefficients kj and 
k 2 , simulated time-response results were collected for this representation. As anticipated, 
due to low dynamic loads and relatively high transmission stiffness, dynamic response 
behaved almost identically to the response of the ideal transmission model with friction 
shown in figures 3.14, 3.15, and 3.16. Consequently, unless high inertial or applied loads 
are experienced by the harmonic-drive system, the inclusion of a compliant element in the 
transmission model provides no significant advantages over the harmonic-drive 
representation with friction alone. However, in harmonic-drive applications designed to 
apply forces or experience high acceleration, the compliant behavior of the transmission 
cannot be ignored and this model may be useful for predicting system behavior in limited 
ranges of operation. Unfortunately, due to kinematic position error, substantial torque 
fluctuations are far too common during normal harmonic drive operation. Because of this 
menacing transmission property, compliant behavior becomes a major factor in determining 

dynamic response and must be addressed for a wide range of harmonic-drive operating 
conditions. 
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3.6 Harmonic-Drive Model 4: 

Transmission with Friction, Compliance 
and Kinematic-Error 

As concluded from the performance of the harmonic-drive model with friction 
losses, prediction of the velocity fluctuations and resonant behavior in the experimental 
response is essential to secure model accuracy in all operating regions. Since velocity 
fluctuations result from kinematic position error, this transmission property must be 
recognized in an improved representation. Also, since resonance occurs when the natural 
frequency of the system, as governed by the inertias and transmission compliance, 
coincides with the frequency of kinematic-error disturbances, the interaction between 
flexibility and position error must also be captured in the model. Lastly, since the 
frequency of kinematic-error fluctuations is directly related to the rotational velocity of the 
transmission, accurate frictional losses in the transmission are required to regulate operating 
speeds. As presented in this section, a new model which incorporates these three 
transmission properties demonstrates improved accuracy over a wider range of system 
operation. 


3.6.1 Description of Model 4 

As illustrated in figure 3.19, the improved harmonic-drive model contains a cyclic 
kinematic-error element as well as the transmission compliance and frictional losses 
discussed in the previous models. Since the kinematic error in the harmonic-drives is 
intimately related to the gear-meshing mechanism itself, it is situated in the model adjacent 
to the gear-reduction element. Although I chose to place this cyclic position-error 
component on the input side of the transmission, it can be reflected to any port of the 
transmission without disturbing simulation performance. This additional component 
introduces a cyclic position function into the model which varies as a function of the 
transmission rotational position. Since kinematic error results primarily from gear-tooth 
placement inaccuracies on the flexspline and circular spline, a strictly accurate 
representation of this error would include two separate error functions, one at the gear¬ 
meshing frequency of the flexspline and the other at the gear-meshing frequency of the 
circular spline. Since these two disturbance frequencies are almost identical, complexity 
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can be reduced with little penalty on model accuracy by combining these two components 
into a single function. 


As discussed in section 2.3.2.4 of this document, a separate cyclic-torque 
disturbance mechanism exists in harmonic-drives which is independent of torque 
fluctuations caused by kinematic error. This phenomenon is caused by radial loads exerted 
by the flexspline on the wave-generator which tend to orient the wave-generator along a 
preferred axis of alignment. Since this input diametral torque-variation can be observed 
when the flexspline is not assembled into the circular spline, it has no relationship to the 
kinematic inaccuracies in gear-tooth meshing. In order to represent this phenomenon, a 
cyclic input torque acting at the wave-generator was included in the model to capture this 
spring-like behavior. However, as simulation results indicated, since the magnitude of this 
torque function was small compared to the kinematic-error disturbance, its affect on 
dynamic performance could only be observed at slow operating speeds if at all. Based on 
this conclusion, this cyclic input torque will be ignored in the model derived in this section, 
but should not be forgotten in future models which require enhanced accuracy. 
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3.6.2 Derivation of Dynamic Equations 

Using the notation outlined in figure 3.20, the constitutive relationships for the 
transmission components that have been described previously can be summarized as 
follows. First, as dictated by the position, velocity and torque constraints specified in 
equations (3.15), (3.16), and (3.17) the ideal three-port gear-reduction element can be 
characterized by three relationships: 


0n_wg - (N + 1) 0 n _cs ' N 0 n _f s , 

(3.28) 

®n_wg — (N + 1) C0n_ cs - N co n f s , and 

(3.29) 

T — 1 T —It 

- wg “ (N + 1) - s ~ N n - fs * 

(3.30) 


Second, using the non-linear damping function developed in model 2, the constitutive 
relationship for the friction element in the model is 


^b ^b_constanl "b Tb_dynamic ~b Tb_cyc]ic > 


(3.31) 
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where the three independent friction-torque functions are defined in equations (3.22), 
(3.23) and (3.24). Third, as discussed in the third model and defined in equation (3.26), 
the torque function provided by the stiffness element can be described by 

Tk = ki (9 W g - 0k_om) + k 2 (0 W g - ©k_om) 3 • (3.32) 

Given these definitions for the gear-reduction, frictional, and stiffness properties of the 
transmission, the constitutive constraints in the model will be complete once the torque and 
rotation behavior of the kinematic-error element is characterized. 

As observed in experimental measurements, kinematic position-error signatures in 
harmonic drives can occur at many different frequencies. However, as observed in 
dynamic response data and verified by model performance, the only frequencies which 
influenced the dynamic behavior of the three different harmonic-drive systems occurred at 
one, two, and four cycles every wave-generator revolution. The reason for this is twofold: 
(1) the magnitude of kinematic error at other frequencies was smaller than these primary 
error components, and (2) other frequencies were either too high or too low to excite 
system resonance for the observed range of harmonic-drive operating velocities. For other 
systems which operate at slower velocities and have lower natural frequencies, for 
example, kinematic error at higher harmonics may need to be included. Based on this 
understanding, the representation chosen to describe the position error in the harmonic- 
drive transmission is a harmonic series of the three error components: 

Qerfn — Aj sin (9 wg + <}>i) + A 2 sin (29 W g + <J> 2 ) + 

A 3 sin (49 wg + <J) 3 ) . (3.33) 

Notice that, in this equation, each term is a sinusoidal function of the wave-generator 
rotation, 9 wg , with a specific amplitude, Aj, and phase, <j>i, describing the magnitude and 
relative orientation of each error component. Since these error components result from gear 
error on the flexspline and circular spline, a strictly accurate representation would have 
these functions vary with harmonics of the gear-tooth meshing frequencies. For example, 
in joint configuration 1 where the flexspline is attached to ground, every rotation of the 
wave-generator propagates the tooth-engagement zone completely around the flexspline 
circumference but fails to reach two teeth on the circular spline. Consequently, the 
frequency of the flexspline gear-tooth error components correspond exactly with the 
rotational frequency of the wave-generator while the frequency of the circular-spline gear- 
error components vary with the difference between the wave-generator and circular spline 
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rotation. However, since the input rotation is much larger that the output rotation in the 
transmission, this distinction can be neglected and the aggregate kinematic error can be 
varied as a function of the wave-generator rotation alone. In models which desire to 
capture the kinematic behavior of the harmonic drive to a higher level of detail, this 
clarification may be important. 

Given this representation for the position error in harmonic-drives, the constitutive 
relationships for the kinematic-error element in figure 3.20 can be defined. First, using the 
error-function, 0 er fh> given in equation (3.33) the position constraint on this element can be 
specified as 

®n_wg = 0k_out + N 0 e rfn • (3.34) 

Since the magnitude of the error function is measured in terms of degrees of output 
rotation, it must be scaled by the transmission gear-ratio. The velocity constraint on this 
element can be defined by taking the derivative of the position constraint: 

“n_wg = W k _cm, + N £ e er f„ , (3.35) 

where, given the relationship in (3.33) the derivative of the error function can be defined as 

cj~ ®erfn = &Hvg Aj COS (0 W g + (j)]) + 2 C0 W g A 2 COS (20 W g + ({> 2 ) + 


4 co wg A 3 cos (40 wg + (j> 3 ), or (3.36) 

^ 0 e rfn = 0 ) wg 0 ) erfh , where ( 3 . 37 ) 

^erfn = Aj COS (0 W g + (j>i) + 2 A 2 COS (20 W g + (^ 2 ) + 

4 A 3 cos (40 wg + (j> 3 ). (3.38) 

Using this definition, the velocity constraint in equation (3.35) can now be rewritten as 

^n_wg = + N 0) W g 0) er f n . (3.39) 


By applying the law of power conservation, this velocity constraint can be used to derive 
the torque constraint across the position-error element. First, the power at the input port of 
the element must always equal the power at the output: 

Pin = Pout > Of 4 Q) 

T k _out ^k_out ~ fn_wg ^n_wg • (3.41) 
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By substituting the velocity constraint in (3.39) into this equation, the torque constraint 
becomes 


Tk_out COk _out ~ T n _ W g (^k_out ■*" N co W g co er f n ). ( 3 . 42 ) 

Since the velocity variations across the spring will be small compared to average velocity 
values, C0k_out can be assumed to be identical to co wg and the torque relationship reduces to: 

Tk_out = T n wg (1 + N COerfn) • (3.43) 

Therefore, using this relationship and the equation for G) erfn in (3.38) the torque constraint 
required to ensure power conservation across the kinematic-error element is complete. 

Having developed the position, velocity, and torque constraints for the position 
error in the model, the constitutive relationships governing model performance are 
complete, and the equations of motion can be finished by specifying the remaining 
kinematic and torque-balance relationships. Specifically, as seen in equations (3.18) and 

(3.19), a torque-balance across the output ports of the harmonic-drive yields the 
relationships 


Tfs = T n _f s - T b , and 
T cs = T n cs - T b . 

Also, by performing a similar force balance across the spring on the input of the 
transmission the following relationship results: 

T w g = T k = T k out . (3.46) 

Lastly, by specifying the kinematic constraints at the output ports of the harmonic-drive, 
namely 


(3.44) 

(3.45) 


0n cs — @cs » and 


0n_fs — 0fs > 


(3.47) 

(3.48) 


the equations of motion are complete. Therefore, using the equations presented in this 
section, the positions and velocities at the three harmonic-drive ports can be directly related 
to the torques across the harmonic-drive. 
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3.6.3 Calculation of Input Parameters 

In order to simulate dynamic performance using the equations for model 4, 
estimates need to be made for the values of the friction, stiffness and kinematic error in the 
harmonic-drives. Using the experimental friction data provided in section 2.5 of this 
document, the constant, dynamic, and cyclic friction values required by equation (3.31) 
were easily determined. Since these values were extracted directly from dynamic response 
observations, they provided accurate and reliable results in the dynamic model. 

Estimating the values for the non-linear stiffness element in model 4 proved to be 
more complex than capturing friction behavior. In particular, due to the precarious 
influence of friction on the static stiffness measurements outlined in section 2.3 of this 
document, the resulting stiffness values delivered poor dynamic performance when inserted 
into the model. As a more direct method of determining the dynamic stiffness properties of 
each transmission, the experimental dynamic-response measurements outlined in section 
2.6 were used to identify the resonant frequencies of the three harmonic-drive systems. 
From this known frequency and the given inertia values, the linear component of the ideal 
stiffness element could be adjusted to match the observed resonance. Given this linear 
component, the cubic stiffness coefficient required by equation (3.32) was selected so that 
the static stiffness profile generated by the model most closely resembled the experimental 
static stiffness behavior. Figures 3.21, 3.22, and 3.23 illustrate the resulting stiffness 
profiles used in the three joint-models as compared to experimental data and catalog 
predictions. In each of these plots, it can be seen that the theoretical stiffness profiles show 
significant deviation for the static stiffness measurements. In particular, the hysteresis loss 
in each curve, as dictated by the constant friction component used in the model, Tb constant, 
while producing good dynamic response, is unreliable as a predictor for static hysteresis 
loss. Additionally, the shape of the theoretical stiffness profiles required to ensure good 
dynamic response in the three harmonic-drive models provides a poor estimate of the 
measured static-stiffness profiles. As this comparison illustrates, for all three joints, the 
experimental measurements typically yield static stiffness values that are higher than the 
actual dynamic stiffness in the drive. This discrepancy is not surprising since Coulomb¬ 
like friction in the transmission can increase the measured torque values during static 
stiffness measurement. To address this problem, a new model will be developed later 
which predicts these Coulomb-friction effects. 
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Figure 3.21: Joint 1 stiffness profile for harmonic-drive model 4 
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Figure 3.22: Joint 2 stiffness profile for harmonic-drive model 4 
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Output Torque (NTm) 


Figure 3.23: Joint 3 stiffness profile for harmonic-drive model 4 


Given the friction and stiffness properties, the only remaining parameters that need 
to be determined are the magnitudes of the kinematic-error components in equation (3.33). 
The amplitudes of the three error components were determined by matching the amplitude 
of the dynamic position error predicted by the model to experimental values. In theory, if 
kinematic error on the testing apparatus could be measured under truly static conditions, 
direct estimates of these error components could be made. However, since data on the 
apparatus can only be collected for non-zero operating velocities, dynamic vibrations which 
induce transmission compliance contaminate the kinematic-error data. As presented in 
section 2.2 of this document, by operating the harmonic-drives at slow velocities where 
resonance vibration was minimal, measurements were made of the kinematic error in the 
three harmonic drives. In table 3.4, these experimental estimates are compared to the error 
values which provided optimal dynamic performance in the theoretical model. By 
comparing these results, it can be seen that, in most cases, values near the experimental 
measurements provided accurate dynamic response simulations, but, in a few instances, 
differences between experimental and theoretical values could deviate by up to a factor of 
two. From this discrepancy, two possible conclusions can be drawn: (1) slow operating 
velocities can significantly influence the magnitude of the position-error measured in the 
experimental apparatus, and (2) the model developed in this section does not fully capture 
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Table 3.4: Kinematic-erro 

r components in 

model 4 vs. measured values 
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the influence of kinematic error on dynamic vibrations in the transmission. Since additional 
dynamic simulations with the theoretical model illustrated that dynamic vibrations could 
cause significant amplification of transmission position error at small operating velocities, it 
is likely that the first of these two conclusions is most responsible for the difference 
between experimental and theoretical parameters. 


3.6.4 Evaluation of Simulated Results 

The harmonic-drive model developed in this section was inserted into the overall 
system models for the two joint configurations and dynamic step-response simulations 
were then performed to replicate the experimental results. Using the experimental 
parameters for friction, stiffness, and kinematic error discussed previously, time-response 
data was collected over a range of motor-current steps for the three harmonic-drive testing- 
station models. These simulated results are presented in appendix E in the identical format 
used for the experimental time-response results of appendix D. In particular, graphs are 
shown which illustrate the position, velocity, and torque response for the observed range 
of harmonic-drive operating behavior. Table 3.5 is reproduced from a similar index in 
appendix E and summarizes the figure numbers for all of the simulated dynamic-response 
plots. By directly comparing these plots to the experimental results in appendix D, the 
capacity of the dynamic model to predict actual harmonic-drive behavior can be evaluated. 

Based on how accurately the theoretical predictions in appendix E compare to the 
experimental results in appendix D, the benefits of this model are summarized below: 

[Tj Position, velocity and torque response in the model shows good agreement to 
actual data at operating speeds removed from regions of resonant behavior. 
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Table 3.5: Summary of 

simulated time-response plots in 

appendix E 
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System resonances occur at the same operating velocities and consequently 
display the same vibration frequencies as experimental results. 

Amplitudes of input and output velocity fluctuation near resonance are 
relatively close to the experimentally observed amplitudes. 


Conversely, by noting the differences between predicted and actual behavior, the 
limitations of this model can also be observed: 





At current-steps which result in velocities that excite system resonance, model 

predictions of operating velocity are typically much higher than observed in 
the experimental results. 

The model also fails to capture the sudden velocity-jump behavior that occurs 
when the system breaks through a resonant region. 

Amplitude of output torque fluctuation near resonance, especially on joints 1 
and 3, is somewhat lower than the experimental torque-fluctuation amplitude. 


These limitations on the accuracy of theoretical predictions lend insight into possible areas 
for model improvement. First, since the predicted velocities near system resonance are 
much higher that the observed velocities, it is clear that this model fails to capture the 
enhanced frictional dissipation which appears due to resonant vibration. As discussed in 
light of experimental results, this dissipation probably results from a Coulomb-like 
frictional torque which increases with the amplitude of torque fluctuations in the 
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transmission. Since, in the three drives tested, the resonance phenomenon dominates 
transmission behavior in a wide portion of the observed operating range* an accurate 
harmonic-drive model cannot ignore these resonance losses. As a second consequence of 
this resonance-dissipation absence, the velocity jump that occurs in experimental response 
when resonance vibrations cease does not appear in the theoretical response. 
Consequently, any improved representation which is able to replicate friction loss at 
resonance is likely to show this unusual velocity behavior as well. Furthermore, since a 
friction model at resonance is likely to supplement dynamic torque fluctuations with 
variations in frictional torque, the low torque-amplitudes observed on the joints 1 and 3 in 
this model may also improve. As a final note, the experimental current response on all 
three joints displays electrical noise which is not reproduced in the current response of the 
model. Since this noise had no effect on the dynamic performance of the actual system, it 
was completely ignored in the simulation. 


3.6.5 Model 4 Conclusions 

The simulated results of the harmonic-drive model presented in this section 
demonstrate that, by including the combined effect of friction, compliance, and kinematic 
error in the model, performance can be significantly improved over models that do not 
acknowledge all three of these transmission properties. In particular, the velocity and 
torque fluctuations which appear in the experimental response can now be reproduced due 
to the kinematic error in the model. Additionally, due to the interaction of the kinematic 
error and compliance, system resonances appear in the theoretical predictions that match the 
experimental results. Unfortunately, since the model fails to capture frictional losses which 
occur at resonance, simulation accuracy decays at velocities which excite resonance 
vibration. Based on these conclusions, it seems likely that this model may be useful for 
predicting performance of dynamic systems that typically operate far from regions of 
resonant behavior. Since the kinematic-error effects are featured in this model, static 
positioning applications which require high accuracy may also benefit from predictions 
delivered by this model. However, for any application which requires accurate predictions 
of transmission behavior over a wide operating range, an improved harmonic-drive 
representation which can predict frictional losses at resonance is required. 

Another conclusion which can be made from this model is that reliable 
measurements of transmission stiffness are difficult to obtain. In particular, when stiffness 
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curves collected from static measurements were used in the model, the resulting dynamic 
response held few similarities to actual behavior. As explained above, this discrepancy can 
be blamed on the inaccuracy of static stiffness measurements due to friction in the 
harmonic-drive. From the substantial differences between the theoretical and actual static 
stiffness profiles discussed previously, it can be concluded that interaction between 
stiffness and friction in the actual system is far more complex than the model recognizes. 
Based on this result, an improved harmonic-drive model, in addition to improving 
predictions at dynamic resonance, should also better represent the static stiffness behavior 
in the transmission. In order to improve model performance, dynamic response data had to 
be used to isolate the natural frequency of the system and derive the corresponding 
transmission flexibility. Although this new stiffness data provided good dynamic results, 
slight changes in its value could cause significant shifts in the operating velocities which 
excited resonant behavior. This unfortunate result demonstrates that, unless an extremely 
accurate model is developed which can reliably predict transmission stiffness from static 
measurements, it is likely that static stiffness data will never be useful for accurate 
predictions of harmonic-drive dynamic behavior. 

Perhaps the most important conclusion that can be drawn from the results in this 
section is that accurate harmonic-drive modeling is a very complex task. The scope of the 
model presented in this section represents the maximum useful complexity that can be 
included in a model that ignores more complex geometric interaction of the harmonic-drive 
components. Since the resulting model performance was completely unreliable for large 
portions of the harmonic-drive operating ranges, the forecast is gloomy for designers who 
wish to make accurate predictions using a simple transmission representation. In a final 
attempt to create a model that can make accurate predictions while still containing a 
manageable amount of complexity, a representation will be developed below which 
incorporates the three transmission properties used in model 4 into a simplified 
representation of harmonic-drive gear-tooth interaction. 
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3.7 Harmonic-Drive Model 5: 

Transmission Including Gear-Tooth 
Geometry 

An accurate representation of harmonic-drive behavior cannot exist without 
incorporating frictional losses, transmission compliance, and kinematic error. However, as 
model 4 illustrated, these three properties are not the only factors which influence 
harmonic-drive behavior. Additionally, a mechanism which dissipates energy during 
system resonance and influences static friction measurements must also be captured. To 
accomplish this task, I have transplanted the friction, compliance, and kinematic error 
effects developed in the previous models into a new framework which imitates the gear- 
tooth meshing action inside the harmonic drive. As discussed in this section, the increase 

in complexity introduced by this new representation is rewarded by a substantial increase in 
model performance. 


3.7.1 Description of Modei 5 

In order to develop a model that mimics gear-tooth meshing action in the harmonic 
drive, a more intimate understanding of harmonic drive operation is required. Figure 3.24 
depicts close-ups of the gear tooth meshing geometry for a section of a harmonic-drive. 
For the case where the flexspline is fixed to ground, clockwise rotation of the wave- 
generator will cause the flexspline gear teeth to be pushed into the circular spline gear-teeth 
along the major axis of the wave-generator. Since the flexspline has two fewer teeth than 
the circular spline, the mating gear teeth will be slightly misaligned before meshing. 
Consequently, as the flexspline teeth are pushed radially into the circular spline teeth, the 
leading, or right-hand-side, face of the flexspline teeth will slide along the trailing, or left- 
hand-side, face of the circular-spline teeth during engagement. The forces generated as the 
teeth slide against each other will cause the gear teeth in the meshing zone to align and, 
consequently, push the circular spline in the clockwise direction. For the case where the 
circular spline is Fixed to ground, clockwise rotation of the wave-generator still pushes the 
flexspline teeth into the circular-spline teeth along its major axis. However, since the 
circular spline cannot move, the forces generated as the flexspline teeth slide along the 
circular spline teeth result in counterclockwise rotation of the flexspline. 
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Using this understanding of gear-tooth meshing action, assumptions can be made 
which reduce the complexity required in a harmonic-drive model which mimics this 
behavior. First, by realizing that during normal harmonic-drive operation, anywhere from 
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ten to greater than fifty teeth can be in contact between the circular spline and flexspline, the 
effect of individual gear teeth can be ignored in a dynamic model. This assumption was 
confirmed by the kinematic error results in section 2.2.2.1.4 which discovered that the 
error of individual teeth was not detectable. Consequently, in order to capture the 
aggregate effects of gear-tooth meshing, a model which describes sliding on a single 
continuous gear-tooth surface should suffice. Second, since the gear-teeth are very small 
compared to the pitch-diameter of the transmission, the localized motion of the gear-teeth 
during meshing shows little trace of the overall rotational motion. Based on this 
assumption, a model which describes this gear-tooth action can be simplified by using 
translational rather than rotational motion. Using these two assumptions, a dynamic model 

that captures the aggregate effect of the gear-tooth meshing mechanism is illustrated in 
figure 3.25. 


The simplified gear-tooth meshing model in figure 3.25 consists of a series of 
massless inclined planes which embody the wave-generator and gear-tooth geometry of the 



Figure 3.25: Schematic of harmonic-drive model with gear-tooth geometry 
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actual transmission. These planes are confined to horizontal and vertical motion by a series 
of ideal rollers and rigid surfaces. In particular, note that the output pons of the wave- 
generator, flexspline, and circular spline are all restricted solely to horizontal motion so 
they can be interfaced easily with an overall system model. Although this cannot be 
illustrated in the schematic, by assuming that all of the surfaces which constrain motion are 
infinitely long, this model can be used to describe all ranges of transmission operation. 

By comparing the actual transmission in figure 3.24 to the model in 3.25, several 
parallels can be drawn between actual and modeled behavior. First, horizontal motion in 
the model corresponds directly to rotation of the harmonic-drive components and vertical 
motion in the model corresponds to radial movement in the harmonic drive. Second, the 
low-angle plane and inclined bearing surface on the bottom of the model schematic 
correspond to the wave-generator in the actual drive. The mechanical advantage and 
kinematic properties of the wave-generator eccentricity are captured by the low-angle of this 
wedge. Since the wave-generator motion pushes the flexspline gear teeth into the circular 
spline, a third parallel can be drawn between the two central wedges and the flexspline in 
the actual system. Similarly, the wedge at the top of the model schematic which mates with 
the flexspline wedge corresponds directly to the circular-spline in the actual system. So 
that this model effectively replicates the gear-tooth rubbing action in the actual 
transmission, the angle of both the flexspline and circular spline mating wedges is identical 
to the traditional harmonic-drive gear-tooth angle of thirty degrees. The angle of the wave- 
generator wedge at the bottom of the schematic can be calculated from the gear-ratio of the 
particular harmonic drive as discussed later. 

A better understanding of the likeness between the model in figure 3.25 and the 
actual transmission can be gained by considering the case where the flexspline is fixed to 
ground. Under this condition, as discussed previously, clockwise rotation of the wave- 
generator in an actual transmission produces clockwise rotation of the circular spline. This 
relative motion results since the wave-generator pushes the flexspline gear teeth into the 
circular spline gear teeth, and the forces at the gear-tooth interface drive the circular spline 
in the clockwise direction. Identical behavior can be seen in the schematic in figure 3.25 if 
the flexspline port is fixed to ground. Specifically, the clockwise rotation of the wave- 
generator can be mimicked by sliding the wave-generator wedge in the model horizontally 
to the right. Due to this motion, the two flexspline wedges in the model, which are not 
allowed to move horizontally, slide vertically upward and push against the circular-spline 
wedge. This vertical motion directly parallels the engagement of the gear teeth that occurs 
in the actual transmission. Continuing the analogy, the vertical motion of the flexspline 
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wedges in the model causes the circular-spline wedge to move horizontally to the right. 
Again, this horizontal motion mimics the clockwise rotation that results in the circular 
spline of the actual transmission. By fixing the circular-spline port in this model to ground, 
a similar exercise can be performed to relate model behavior to the meshing action in an 
actual harmonic drive under a different operating configuration. 

Given this geometric framework which describes the motion and force constraints 
present in harmonic drives, remaining dynamic effects can be incorporated by introducing 
friction, compliance, and kinematic-error properties. Specifically, since it is likely that 
most frictional loss in harmonic-drives is generated from rubbing at the gear-tooth 
interface, all friction in this model is localized at the contact surface between the flexspline 
and circular spline wedges. However, unlike the previous model, Coulomb friction 
properties are also included at this contact surface to provide a mechanism for 
demonstrating increased losses at system resonance. Compliant behavior, on the other 
hand, is not as easy to pinpoint as friction. In particular, past research has indicated that, in 
addition to deflection in both the flexspline and circular spline gear-teeth, significant 
deflection may also occur radially in the wave-generator as well. Since no information was 
available to supervise the distribution of the overall transmission compliance into these 
individual locations, a single stiffness element was included at the flexspline for simplicity. 
In order to maximize the influence of friction at the gear-tooth interface on static stiffness 
measurements in the model, this spring was oriented vertically rather than horizontally. 
Lastly since kinematic error in harmonic drives is caused mainly by gear-tooth 
manufacturing errors on both the circular spline and flexspline, two separate error function 
should be included which vary the positions of the flexspline and circular-spline gear-tooth 
surfaces. However, since the combined effect of these two error elements can be closely 
approximated by a single error function, only one kinematic-error element is included in the 
model. For the sake of simplicity, this element is located at the flexspline. Also, since 
numerical problems can arise if the velocity at the input of this element becomes zero, the 
kinematic error is oriented vertically rather than horizontally . If future models require 

improved accuracy, this error function can be distributed into different locations in the 
model. 

As illustrated in experimental dynamic-response results, such as figures D.1.7, 
D.2.7, and D.3.7 in Appendix D, torque variations in harmonic-drives typically fluctuate to 
large positive and negative values. From the model schematic in figure 3.25, it can be seen 
that, if the normal force at the gear-tooth contact surface becomes negative due to this 
alternating torque, the coulomb friction function at that surface becomes undefined. By 


144 



Chapter 3: Harmonic-Driv e Model 5 


better understanding how the load is carried by the harmonic-drive gear-teeth during these 
torque oscillations, a new model can be developed to capture this behavior. In particular, 
as figure 3.26 illustrates, when the driving torque across the gear-teeth is positive, the load 
is carried by the gear-teeth surfaces which push the circular spline in the positive direction. 
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However, when the driving torque becomes negative due to inertial vibrations, the 
transmission velocity is still in the same direction but the opposite face of the gear teeth 
must bear the resulting load. In order to capture this behavior in a dynamic model, the 
representation shown in figure 3.27 can be used whenever the normal force at the gear 
tooth interface becomes negative. Since this new model is essentially a mirror image of the 
previous one, the normal force at the gear-tooth interface becomes positive and resulting 
coulomb friction can be easily determined. Using this representation along with the 
previous one, an attempt was made to mimic both positive and negative gear-tooth loading 
in the harmonic-dnve model by rapidly switching between models during simulation. 
Unfortunately, since the numerical simulation used to solve the equations of motion was 
not very tolerant of discontinuous changes in the dynamic equations, a solution could not 
be reached. Instead, the single representation in figure 3.25 was used exclusively, and 
whenever the normal force at the gear-tooth surface became negative, the coulomb friction 



Figure 3.27: Revised transmission model for negative-torque loading 
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was set to zero. This inaccurate assumption most likely represents the main shortcoming of 
this model and identifies an important area for further improvement. 


3.7.2 Derivation of Dynamic Equations 

Before the derivation of the equations of motion for this system can begin, notation 
conventions must first be established. Specifically, the variable names for the positions 
and velocities and their corresponding sign conventions are illustrated in figure 3.28. As 
this figure illustrates, positions are defined by the symbol 0 and velocities are represented 
by the symbol 0 ). Sign conventions for all of the kinematic variables in the model, except 
0f s and G)f s , are defined to be positive in the direction they will travel when the wave- 
generator is moved in the positive direction. Similarly, torque conventions are illustrated in 
figure 3.29. For all of these variables the symbol T is used, and positive torque is defined 
in the direction it will act when a load is applied in the positive direction on the wave- 
generator while the flexspline and circular spine are kept stationary. As illustrated in both 
these figures, the geometric angles of both the wave-generator and the gear-teeth in the 
model are defined by the symbol a. These geometric angles should not be confused with 
the angular position coordinate, 0. 


Given the variable names provided in these two figures, the geometric constraints in 
the model will be used to define the kinematic equations governing model motion. First, 
due to the inclined face of the wave-generator wedge, the horizontal motion of the wave- 
generator and flexspline can be directly related to the vertical motion of the flexspline base 
by the equations 


®tooth_base tan otwg (®wg ~ 0fs) •> and 


(3.49) 


t^tooth_base tan Ot^g (CDwg - COfs) • 


(3.50) 


Similarly, due to the gear-tooth angle, the circular-spline and flexspline horizontal motion 
can be used to determine the vertical motion of the flexspline tip: 


®tooth_tip v^cs ” 0fs) > &nd 

tan cXtooth 

(3.51) 

^tooth_tip (tOcs C0f s ) . 

tan (Xt 00 th 

(3.52) 
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Using these equations and the geometry of the model, the equations specifying the motion 
at both the wave-generator surface and the gear-tooth surface can be derived: 

0\vg_surface = COS (X wg 0 wg + sin 0C wg 0 tO oth_base ' COS (X wg 0f s , (3.53) 

W wg_surface = COS 0C wg C0 wg + sin a wg CDtoo^basg - COS 0t wg G)f s , (3.54) 

0tooth_surface = COS CC t00t h 0 tO oih_tip + sin a toot h 0 CS - (3.55) 

sin 0 Ct oo th 0f s , and 

w tooth_surface = COS OC toot h CO tooL h_ t i p + sin OC tool h C0 cs - (3.56) 

sin ottooth C3f s . 


Now that the kinematic relationships are complete, the angle of the wave-generator wedge 
required to reproduce the motion in the actual transmission can be found. Specifically, 
given the gear-tooth angle, Ottooth> and the catalog gear ratio, N, the wave-generator angle 
can be determined exactly from the equation: 


a wg = arctan 


(N + 1) tan a t00 ih 


(3.57) 


If these kinematic constraints are indeed accurate, they should produce the same 
relationship between the positions and velocities of the wave-generator, flexspline, and 
circular spline as observed in the actual system. As a confirmation of this fact, while 
ignoring kinematic error and compliance in the model, equations (3.49), (3.50), (3.51), 
(3.52), and (3.57) can be combined to reproduce the ideal three-port transmission 
equations presented in section 3.3. Based on this result, the position and velocity 
relationships presented above provide a complete and accurate description of the kinematic 
behavior of harmonic drives. 


Now that the kinematics of the harmonic-dnve model with gear-tooth geometry are 
better understood, the constitutive relationships which govern the behavior of the individual 
model elements will be derived. First, as seen in previous models, the friction at the gear- 
tooth interface is composed of (1) a velocity-independent, constant torque, (2) a velocity- 
dependent, damping torque, and (3) a cyclic torque function. However, as a new and 
important improvement to this model, a Coulomb friction function is also included at the 
gear tooth surface. The resulting aggregate friction function is 

= Tb_constant + Tb_dynamic + Tbcyclic + Tb_ CO ulomb , (3.58) 
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where 


Tb_constant = ^constant , 

(3.59) 

Tb_dynamic — ^1 ^tooth_surface ^2 (^tooth_surface)^ > 

(3.60) 

Tb-cyclic “ Ab ((®cs “ 0fs) ^b) > ^nd 

(3.61) 

Tb_coulomb ~ M- Ttooth_normal • 

(3.62) 


In these equations, p. represents the Coulomb friction coefficient at the gear-tooth interface 
and the remaining parameters, as seen before, are defined in table 3 . 3 . 

Unlike this friction relationship which contains a new term for the Coulomb losses, 
the stiffness equation used in this model resembles closely the cubic relationship seen in 
previous models: 

Tk ~ ki (0err " 0tooth_tip) ^2 (0err “ ®tooth_tip) • (3.63) 

Following the example set by model 4, the kinematic error function used in model 5 
will be a harmonic series of three error components. However, slightly different than 
model 4, since the kinematic-error function is located on the flexspline, it should vary as a 
function of the flexspline gear-tooth meshing position rather than the wave-generator 
position. Considering that every tooth on the flexspline will mesh once when the wave- 
generator makes one complete revolution relative to the flexspline, the position which 
governs the kinematic-error function should be 0 wg - 0 fs . Adopting the syntax used for 
model 4 in equation (3.33), this kinematic error function becomes 

Gerfn = Ai sin (0 wg _rei + <h) + A 2 sin (20 wg _ re i + 62 ) + 

A 3 sin (40 wg _ rel + 4> 3 ) , (3.64) 

where 

®wg_rel = 0wg ’ 0f s • (3.65) 

Similarly, the derivative of this function can also be defined as 

^ 0erfn = co wg _ re i CD er fh > where (3.66) 
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w erfh — Ai COS (9 W g_rel + 0]) +2 A 2 COS (20 wg _ re i + ({> 2 ) + 

4 A 3 cos (40 wg rel + 4 ) 3 ), and (3.67) 

^wg_rel = W wg - COf s . (3.68) 

Note that, unlike previous models, since the magnitude of the kinematic-error function in 
model 5 does not translate directly to rotation of the flexspline or circular spline, the values 
of the position-error amplitudes, Aj, A 2 , and A 3 , must be adjusted as discussed later. 
Given the position-error function and its derivative, the position, velocity and torque 
constraints can be derived as done for model 4. Specifically, by definition, the position 
and velocity constraints across the kinematic-error element are 


9err ~ 9tooth_base + Qerfn > & n d (3.69) 

w err = ^tooth_base + &>wg_rel &Wn • (3.70) 

By requiring that power be conserved, the following equation can also be derived: 

Terrjn w tooth_base = Tk (0 err . (3.71) 

Using the kinematic-error relationship in equation (3.70) and the geometric constraint in 
equation (3.50), this power relationship can be converted easily into the desired torque 
constraint on the kinematic error element: 


T 


err in ~ 


1 + _9kfh_j x 
tana wg ; 


(3.72) 


At last, with the torque, velocity, and position relationships for the kinematic error fully 
defined, the constitutive relationships required in the model are complete. 

Having calculated both the kinematic and constitutive constraints on the model, the 
task that remains is to balance the forces across the elements in the transmission. This can 
be accomplished by summing the forces separately for the four different geometric shapes 
in the model. For the wedge representing the wave-generator, balancing the forces in the 
horizontal and vertical directions yields the equations 



Twg.normal sin OC 


wg , 


and 


(3.73) 


Tinput_normal — T 


wg_normal COS a W g 


(3.74) 
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Likewise, the equations for the flexspline base are 

T fs = T fs _ .normal “ Twg_ norma i sin (X W g , and 

Terr_in = ^wg-jiormal COS OC W g , 

and the equations for the flexspline tip are 


(3.75) 

(3.76) 


T f s_normal~ Ttooih_normal COS CXtooth" Tt 0 oth_friction sin OC tool h, and (3.77) 

— Ttooth_normal sin OCt 0 oth 7 , tooth_friction COS Ot^oth’ (3.78) 

Using the same procedure, the force-balance equations for the remaining circular-spline 
wedge complete the dynamic constraints: 

Tcs = T t0 oth_normai cos oc toot h - T toot h_f r i ct i on sin oc toot h , and (3.79) 

7'output_normal — 3\ooth_normal Sin OC tool h+ Ttoolh_friction COS CCtoot^. (3.80) 

To solve these equations, (3.74), (3.76), (3.78), and (3.80) can be combined with the 
constitutive torque-relationship for the kinematic error element in order to determine 
Tinput_normai> 3\vg_normai> ^outpu^normaij und T^ 0 Q^^ n Qj- ma j in terms of known quantities. 
Then, by substituting these values along with equation (3.77) into equations (3.73), 
(3.75), and (3.79), the expressions for the torques on the wave-generator, flexspline, and 
circular spline can be found. When the positions and velocities of the three harmonic-drive 
ports are known, these equations, combined with the kinematic and constitutive 
relationships presented above, can be used to calculate all torques experienced by 
transmission. 


3.7.3 Calculation of Input Parameters 

Before the equations in model 5 can be used to simulate dynamic response, values 
for friction, stiffness and kinematic-error properties need to be determined. In the case of 
friction, experimental measurements presented in section 2.5 provided good results when 
used in previous dynamic models. However, in order to use these same values in model 5, 
the experimental friction results, which were measured from the input side of the 
transmission, need to be reflected to the gear-tooth interface, where they will appear in the 
model. The appropriate scaling factor which relates the wave-generator velocity to the 
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gear-tooth surface velocity can be easily found from the kinematic relationships presented 
above. On joint configuration 1, for example, this scaling factor is 


B = 


_ 1 _ 

(N + 1) sin atooth 


, where 


(3.81) 


^tooth surface — B 03 


wg 


(3.82) 


Given this scaling factor, the values for the tooth-surface friction coefficients required to 
generate the same losses observed experimentally can be calculated readily. By assuming 
that any power dissipated by a friction source located at the wave-generator should be 
identical to the power dissipated when the friction is located at the gear-tooth interface, the 
experimental friction results can be directly related to friction parameters at the gear-tooth 
interface. In the case of velocity independent friction, this analysis yields the result 


t _ u constant_wg 

Constant-^- , 

D 


(3.82) 


where b constant is the constant friction at the gear-tooth interface from equation (3.59), 
bconstant_w g is the constant friction measured at the gear-tooth interface, and B is the scaling 
factor presented in equation (3.81). Since the amplitude of the cyclic friction in equation 
(3.61) is also independent of velocity, it can be calculated in an identical manner. 
However, the dynamic friction coefficients in equation (3.60), since they are linear and 
cubic functions of velocity, must be scaled by factors of B 2 and B 4 , respectively. 

Since no experimental data were collected to measure the Coulomb friction in the 
transmission, the remaining friction equation, (3.62), cannot be calculated using the same 
method used to derive the other friction terms. Instead, tabulated handbook values for 
typical Coulomb friction coefficients were used. Specifically, for lubricated steel-to-steel 
contact, listed values ranged from 0.03 to 0.16, and using this range as a guide, friction 
coefficients were then adjusted to optimize model performance. The resulting coefficients 
that delivered good model performance were 0.15 for joint 1, 0.05 for joint 2 and 0.10 for 
joint 3. However, simulated results illustrated that small changes in these friction values 
within the quoted range could greatly influence the dynamic response. This unfortunate 
result identifies a serious barrier to accurate harmonic-drive modeling since reliable 
Coulomb friction values are notoriously difficult to identify. 


Having determined useful values for the friction in the harmonic-drive model, the 
second property that needs to be evaluated is transmission flexibility. As was the case in 
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the previous model, static stiffness measurement were of little help in determining the ideal 
compliance for the model. Instead, the resonances that appeared in the experimental time- 
response data were used to identify the natural frequencies of each harmonic-drive system. 
Given these frequencies, the inertias seen by the spring, as reflected through the geometry 
of the model, were used to calculate the linear component of the ideal compliance required 
to produce the same resonant behavior. With the linear coefficient fixed, the cubic stiffness 
coefficient was adjusted to reproduce the curvature in the measured static-stiffness profiles 
and provide good dynamic-response behavior around resonance. In particular, cubic¬ 
stiffness values were adjusted to increase or decrease the velocity range over which 
vibrations were sustained in each system resonance. Using these stiffness parameters, 
static stiffness simulations were performed using the harmonic-drive model, and results are 
shown in figures 3.30, 3.31, and 3.32. 

By comparing the theoretical stiffness predictions to experimental results and 
catalog values, as shown in the three figures, useful insights about model integrity can be 
gained. First, since the Coulomb friction introduced in this model influences static 
stiffness predictions, the simulated theoretical compliance curves show a stiffer profile than 
the ideal stiffness alone would predict. More specifically, under static loading conditions, 
in order to effect torsion at the transmission output, both the model stiffness and Coulomb 



Figure 3.30: Joint 1 stiffness profile for harmonic-drive model 5 
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Figure 3.31: Joint 2 stiffness profile for harmonic-drive model 5 
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friction must be overcome. Since the magnitude of the Coulomb friction force at the gear- 
tooth surface increases as the load on the surface increases, the applied load required to 
overcome this friction force increases linearly as the load increases. During static stiffness 
measurements, this variation in Coulomb-friction force will appear as an increase in the 
measured stiffness during joint loading and a decrease in stiffness when the joint is 
unloaded. This behavior can be observed in all of the theoretical stiffness profiles in 
figures 3.30, 3.31, and 3.32. As illustrated most clearly on joints 1 and 3, this apparent 
increase in stiffness during loading produces theoretical stiffness measurements that more 
closely resemble the experimental stiffness profiles. This result is positive evidence that the 
gear-tooth mechanism captured in the model more closely mimics the physical behavior in 
the actual system. However, since in all cases, the overall shapes of the predicted stiffness 
profiles show poor agreement to measured values, the actual friction mechanism in 
harmonic drives is probably much more complicated than this model recognizes. 

The second hallmark which can be used to evaluate the accuracy of the theoretical 
stiffness predictions in hysteresis loss. In the model, hysteresis in the calculated stiffness 
profiles is governed by the velocity-independent friction component as well as the Coulomb 
friction at the gear-tooth interface. In particular, the constant friction term provides a fixed- 
width difference between the loading and unloading stiffness curves while the Coulomb 
friction adds a hysteresis-width component that increases linearly with the applied load. 
The combined effect of these two friction factors can be seen in the theoretical stiffness 
curves in figures 3.30, 3.31, and 3.32. From the poor agreement between the predicted 
and experimental hysteresis loss it can be concluded that, although the specified friction 
values provided good dynamic response, the interaction between friction and compliance in 
the actual system is more complex that the model acknowledges. 

Now that values for the friction and compliance effects in the transmission model 
have been isolated, the only parameters that remain unknown are the kinematic-error 
amplitudes of equation (3.64). As discussed previously, since low-velocity measurements 
of these position-error components proved somewhat unreliable due to dynamic vibrations, 
the values were selected in the previous model by matching predicted and experimental 
dynamic position-error. In theory, since the resulting values used in model 4 should 
describe the absolute kinematic error in the transmission, the same numbers should also 
apply to the new harmonic-drive model with gear-tooth geometry. When these kinematic- 
error amplitudes were reflected appropriately through the geometry of model 5 and 
incorporated into the simulation, however, the resulting dynamic-response predictions fell 
short of performance expectations. Specifically, the torque-fluctuation amplitudes that 
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resulted from the given kinematic-error values were much smaller than the experimental 
values. Because of these diminished torques, the amount of dissipation at resonance due to 
Coulomb friction was lower than observed in the actual system. To remedy this problem, 
it was discovered that, with slight adjustments to the calculated kinematic-error values, 
model predictions improved greatly. The final values which delivered optimal model 
performance are listed in Table 3.6 and compared to the experimental predictions from 
section 2.2. Based on this evolution of reliable kinematic-error values, a few useful 
conclusions can be derived. First, since the amplitudes which provided good results in 
model 4 performed less admirably in model 5, it is clear that the mechanisms used to 
introduce kinematic-error effects in both these models are incomplete. For example, as 
discussed previously, if the model including gear-tooth geometry could also incorporate the 
transmission behavior when the normal force at the gear-tooth interface reverses direction, 
it is likely that the influence of kinematic error on dynamic behavior could change 
dramatically. Additionally, from the difference between low-velocity measurements and 
theoretical values shown in table 3.6, the reliability of kinematic-error measurements which 
are not taken under truly static conditions is further questioned. In conclusion, the 
unfortunate result which can be derived from these observations is that the accurate 
kinematic error estimates are elusive. However, through careful parameter adjustment 
based on a good understanding of system dynamic operation, values can be chosen that 
deliver relatively good model performance. 


Table 3.6: Kinematic-erro 

r components in 

model 5 vs. measured values 
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0.020 

0 008 

iiiiliP. i axperi mental 

0.001 

0.013 

0.007 

4omt 2 : theoretical 

0.003 

0 022 

0.006 

0.004 

ipiiiiill • experimental 

•.’•r*x*:’:»;’;*!’:*j*:*:VA’!»r*ti;W:^yt’;*:*:*:':*;»:*j****‘** , ***’*******'i , J’»'^*-*!'i*i*i*!*!*!’:*:*:*j*:*:*******’**' , » , **i** , *'^'*^**t’!*iv:': 

0.001 
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0.003 
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0.007 
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3.7.4 Evaluation of Simulated Results 


By inserting the harmonic-drive model developed in this section into the overall 
system models for the two joint configurations, dynamic step-response simulations were 
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performed to replicate the experimental results. Using the experimental parameters for 
friction, stiffness, and kinematic error discussed previously, time-response data was 
collected over a range of motor-current steps for the three harmonic-drive testing-station 
models. These simulated results are presented in appendix F in the identical format used 
for the experimental time-response results of appendix D and the model 4 simulation results 
in appendix E. In particular, graphs are shown which illustrate the position, velocity, and 
torque response for the observed range of harmonic-drive operating behavior. Table 3.7 is 
extracted from a similar index in appendix F and summarizes the figure numbers for all of 
the simulated dynamic-response plots. By directly comparing these plots to the 
experimental results in appendix D, the capacity of the dynamic model to predict actual 
harmonic-drive behavior can be evaluated. As an explanatory note, the electrical noise 
present in current-sensor readings was not reproduced in the model simulation since it had 
no influence on the dynamic behavior of the system. 


Table 3.7: Summary of 

simulated time-response plots in 

appendix F 
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From these results, several performance advantages of the harmonic-drive model 
with gear-tooth geometry can be noted: 




Due to Coulomb friction losses at the gear-tooth interface, the steady-state 

velocities reached by current-steps which excite system resonance are much 
closer to the experimental values than in previous models. 

The surprising jumps in velocity that appear in the experimental response 

when the velocity exceeds the resonance range are also observed in the 
simulated response. 
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0 Velocities which excite resonance behavior are almost identical between 
predicted and actual response. 

0 As seen previous models, overall predicted response shows good 

agreement to experimental data at operating speeds removed from regions of 
resonant behavior. ‘ * °* 


By noting the difference between the experimental and theoretical results, the shortcomings 
of this model can be recognized as well: 

0 On joints 2 and 3, in particular, input currents that could generate velocities 

high enough to break through resonant regions were lower than observed in 
experimental response. 

0 Compared to experimental response, amplitudes of torque fluctuation were 
too low on joint 1 and too high on joints 2 and 3. 

0 The predicted dynamic position error at system resonance was greater than 
experimental observations for all three joints. 

From the improvements in simulated response shown by this model over its predecessors, 
it is clear that gear-tooth geometry and Coulomb friction are essential components for 
predicting harmonic-drive behavior during resonance. However, from the discrepancies 
between the theoretical and actual response, it also clear that a more detailed representation 
is required before reliable predictions of harmonic-drive performance can be made. In 
particular, as discussed previously, a model that can accommodate the shift in system 
dynamics that occurs when the torque at the gear-teeth changes direction should be able to 
improve predicted results. Currently, the model in this section avoids this complication by 
setting the Coulomb friction at the gear-tooth interface to zero whenever the normal force 
on the gear-tooth surface becomes negative. It is likely that a new model, which captures 
more accurately the actual friction behavior at the gear-teeth, will improve predicted 
performance by (1) increasing loss at system resonance and (2) limiting dynamic position 
error and velocity fluctuations while preserving torque amplitudes. If the new model can 
indeed deliver these effects, the model predictions should resemble more closely the actual 

dynamic response and, consequently, improve the accuracy and reliability of the 
representation. 


By analyzing the influence of model parameters on overall dynamic response 
behavior, the sensitivity of a few important variables was characterized. First, variations in 
transmission compliance, which can influence greatly the natural frequency of vibration, 
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caused substantial deviations in the resulting dynamic response. Since experimental results 
illustrate that resonant behavior governs system performance for a wide range of operating 
velocities, it is not surprising that small variations in model stiffness can have a profound 
effect on dynamic response. Second, variations in the coefficient of Coulomb friction, \i, 
were also shown to influence dynamic response especially at system resonance. 
Surprisingly, decreasing the value of this coefficient often resulted in increased resonance 
losses due to heightened system vibration and torque fluctuations. Lastly, the output 
friction included in the system model had a surprisingly important relationship on dynamic 
behavior. Because increases in this variable resulted in higher normal-forces at the gear- 
tooth interface, resonance losses from Coulomb friction were very dependent on its value. 
A model in which Coulomb friction acts for both positive and negative normal torques may 
not show this same dependence. These observations emphasize the necessity of accurate 
characterization of harmonic-drive stiffness and friction properties for reliable model 
performance. 


3.7.5 Other Models 

In order to evaluate alternative strategies which might improve the performance of 
the harmonic-drive model with gear-tooth geometry, additional representations were 
developed. First, in order to capture the compliance of the torque sensors on each joint, an 
additional spring was included in the model at the grounded port of the transmission. 
However, because of the kinematic constraints of this new compliance and the cubic 
relationship of the harmonic-drive flexibility, the algebraic equations of motion could not be 
solved analytically. To remedy this problem, an additional state variable was introduced 
into the model by attaching an inertia to the spring on the grounded port of the 
transmission. This additional variable allowed for easy algebraic solution of the resulting 
equations, but, since the system was very stiff, numerical solutions were time-consuming. 
Results collected from a few simulations, however, indicated that the new inertia and 
flexibility on the grounded port of the harmonic-drive had a negligible effect on the 
dynamic response of the system. Finally, as mentioned previously, a numerical strategy 
was developed to solve simultaneously the two different sets of dynamic equations that 
describe harmonic-drive operation for positive and negative torques at the gear-tooth 
interface. Unfortunately, since no numerical algorithm was found that could solve the 
discontinuous set of equations reliably, I will resign this problem for later research. 
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3.7.6 Model 5 Conclusions 

The model presented in this section represents a first attempt to incorporate gear- 
tooth-rubbing action into a dynamic model of a harmonic drive. As the simulated results 
illustrated, since this model included Coulomb friction losses at the gear-tooth interface, 
energy dissipation increased during system resonance and dynamic response integrity 
improved greatly in these areas. However, even with this new level of detail in the model, 
discrepancies between predicted and actual response indicate that the dynamics of 
harmonic-drive operation are still even more complex. However, since this is the only 
model which can boast even partially describing harmonic-drive response in all ranges of 
operation, it may be useful in any application in which understanding harmonic-drive 
operation through resonant regions is required. Since the mathematical complexity of this 
model is relatively minimal, it is an ideal candidate for making quick estimates of operating 
behavior, but, if more accurate information is required, I hope that this model can provide a 
template for more complex representations which better capture the subtleties of harmonic- 
drive behavior. In the very least, the model presented in this section can prescribe a 

definite lower bound on the level of detail required in harmonic-drive models to deliver 
reliable results. 

One of the most discouraging conclusions that can be drawn from the results in this 
section is that accurate characterization of transmission properties can be very difficult but, 
unfortunately, also essential for good model performance. In fact, the only way in which 
accurate values could be determined for the stiffness and coulomb friction in this model 
was through careful measurement of dynamic response over a wide range of operation. By 
comparing these dynamic response measurements to the simulated results, the stiffness and 
friction values could be painstakingly adjusted to optimize model performance. Static 
stiffness tests and catalog predictions provided little guidance in this task. Based on this 
result, it seems unlikely that the model in this section, or any similar representation for that 
matter, will ever be able to make accurate predictions of dynamic performance using solely 
catalog values or simple experimental observations. Consequently, for applications in 
which an accurate harmonic-drive model is required,there are no shorcuts around a detailed 
experimental analysis of the actual harmonic-drive system. 

Perhaps the most useful feature of the model in this section is its intimate conceptual 
relationship to the actual behavior in harmonic drives. Since this model depicts the 
interaction of the three different transmission components, harmonic-drive properties such 
as friction, compliance, and kinematic error can be placed at locations which have a definite 
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physical meaning in relationship to the real system. Furthermore, because of this 
conceptual analogy, the model can be improved easily and expanded to incorporate new 
properties resulting from future investigation into transmission behavior. Because of these 
features, harmonic-drive designers might find this model valuable for collecting quick and 
useful information about how adjustments in operating parameters, such as gear-tooth 
angle and frictional losses, can influence overall dynamic response. Engineers and 
designers, on the other hand, might find this model valuable for fostering an intuitive 
understanding of harmonic-drive operating behavior. 
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As presented in this document, my research was targeted at understanding the 
operating behavior of harmonic-drive transmissions and developing mathematical 
representations to describe the observed performance. From my experimental 
investigation, several properties of harmonic drives where characterized and related to the 
overall dynamic behavior of the transmission. By developing models to simulate the 
experimental observations, the importance of many harmonic-drive properties was 
categorized and the complexity of the harmonic-drive operating mechanism was exposed. 
This section summarizes the important insights gained from both the experimental and 
theoretical investigations and presents recommendations for future harmonic-drive research 
and implementation. 


4.1 Experimental Conclusions 

The first set of experimental tests I performed were directed at understanding the 
magnitude and origins of kinematic position-error in harmonic drives. From the collected 
results, I confirmed past Findings which identified primary error components to vary 
cyclically at twice the speed of the wave-generator and subsequent harmonics. Further 
exploration of kinematic error at higher frequencies demonstrated that the kinematic 
signature of individual gear-teeth was too small to be measured. Lastly, by isolating 
individual error components at both the circular spline and flexspline gear-meshing 
frequency, I confirmed that the dominant source of kinematic inaccuracy in harmonic drives 
can be traced to the gear-tooth errors. From evidence discovered through dynamic 
modeling analysis, I realized that, due to fluctuating inertial forces which induce 
transmission compliance, kinematic-error measurements can be corrupted greatly even at 
small operating velocities. 
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The second experimental voyage I undertook was the accurate characterization of 
harmonic-drive stiffness properties. In addition to cataloging the non-linear stiffness 
profiles for three harmonic-drives, perhaps the most important insight I gained on this 
journey was that stiffness measurements taken under static conditions are highly influenced 
by static friction in the transmission. Because of this unfortunate result, static stiffness 
measurements were of little use in determining stiffness properties for harmonic-drive 
dynamic models. To make matters worse, these stiffness measurements were observed to 
deviate from catalog predictions and vary widely with wave-generator angle and gear-tooth 
preload. Since the accuracy of harmonic-drive representations developed later was highly 
sensitive to stiffness values, the volatile nature of harmonic-drive stiffness identifies a 
primary barrier to harmonic-drive modeling efforts. 

Static and dynamic friction in harmonic drives were the third suspects to be 
experimentally audited. In particular, starting-torque measurements were made on three 
harmonic drives which confirmed past results demonstrating that factors such as wave- 
generator angle and gear-tooth preload can greatly affect static friction. Next, experimental 
tests performed at non-zero velocities illustrated that the dynamic friction in harmonic 
drives has three components: (1) a velocity-independent torque probably due to coulomb¬ 
like friction resulting from gear-tooth preload, (2) a non-linear dynamic torque caused by 
the lubrication effects in the transmission, and (3) an increased energy loss during system 
resonance resulting from high torque fluctuations which increase coulomb-like friction at 
the gear-teeth. The increased losses at resonance, in particular, were sighted for causing 
unusual and unpredictable behavior during harmonic-drive operation. Lastly, to add an 
extra element of complexity, the velocity-independent friction torque was observed to vary 
considerably with the output rotation of the transmission. 

Using the results collected for harmonic-drive kinematic error, compliance, and 
friction properties, the overall dynamic response of the transmission was scrutinized. 
From this investigation, several behavior patterns common to the three harmonic drives 
tested were identified and related to known transmission properties. First, all dynamic 
response showed significant velocity and torque fluctuations due to the kinematic 
inaccuracies of the transmission. Second, due to the relatively low stiffness of harmonic 
drives, these fluctuations excited dynamic resonances over a substantial range of the 
harmonic-drive operating envelope. Torque and velocity oscillations were observed to 
increase dramatically in these resonant regions. Third, since coulomb-like friction 
increased energy dissipation at resonance, large increases in motor current produced little or 
no increase in operating velocity around resonance. Additionally, due to this energy loss at 
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resonance, dramatic jumps in velocity were observed when the system finally exceeded 
these resonant regimes. Fourth, since the harmonic-drive stiffness is non-linear, the 
natural frequency of vibration increased with increased torque levels at resonance. This 
rise in resonant frequency increased the velocity range which excited system resonance and 
accentuated the velocity jump that occurred when these resonance ranges were exceeded. 
Lastly, due to friction variation with output rotation, the location of the velocity jump at 
resonance exit was observed to vary unpredictably. 

If a single conclusion must be derived from the experimental results of my research, 
it should be that harmonic-drive properties are very dependent on operating and 
environmental conditions. In particular, stiffness and friction investigations demonstrated 
that factors such as input orientation, output orientation, gear-tooth preload, and operating 
temperature can greatly influence measured values. In an attempt to unify these results, I 
propose a final description to explain the observed behavior. Specifically, due to the 
slightly different gear-tooth spacing on the circular spline and flexspline, a small amount of 
gear-tooth deflection is required to ensure tight meshing. This deflection, in turn, dictates 
the necessity for a non-zero preloading force at the gear-tooth interface. As a fortunate 
side-effect of this gear-tooth loading, backlash in harmonic-drives is minimal. However, 
since this force increases the amount of Coulomb friction during gear-tooth rubbing, the 
friction in the transmission is greatly influenced by the tightness-of-fit at the gear teeth. 
Similarly, since harmonic-drive stiffness is probably directly proportional to the metal-to- 
metal contact area at the gear-tooth interface and a tight gear-tooth fit results in a large 
contact area, it is reasonable to assume that harmonic-drive stiffness should vary with gear- 
tooth preloading as well. Using this description of the gear tooth mechanism, an 
explanation can be made to describe the increases in both stiffness and friction that were 
observed in the experimental results when the preload was increased. Developing this 
theory further, the influence of kinematic error on gear-tooth preloading can also be related 
to variations in harmonic-drive properties. Specifically, as kinematic-error results 
demonstrate, geometric inaccuracies exist in the gear-teeth of harmonic drive transmissions. 
These inaccuracies are likely to cause variations in the tightness of fit between the circular 
spline and flexspline gear teeth when these two components are rotated relative to each 
other. Based on this assumption, it is likely that both friction and stiffness in the 
transmission should also vary due to this fluctuation in gear-tooth preload that occurs when 
either the input or output port of the transmission is rotated. Using this explanation, the 
experimental variations in stiffness and friction with input and output orientation can be 
justified. Lastly, since other conditions such as gear-tooth wear and operating temperature 
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can influence the tightness of gear-tooth meshing, stiffness and friction are vulnerable to 
these influences as well. In general, using this argument, any variation in an environmental 
or operating conditions which can affect loading at the gear-tooth interface will be 
manifested in the stiffness and friction of the transmission. I am optimistic that future work 
in this area can not only quantify this interdependence of transmission properties but can 
also develop mathematical representations to predict experimental measurements under 
static and dynamic conditions. 


4.2 Theoretical Modeling Conclusions 

In my theoretical investigation, five different harmonic-drive models are presented 
to illustrate the level of detail required to reproduce actual harmonic-drive behavior to 
different degrees of accuracy. Since each of these models boasts three-port behavior, all 
possible harmonic-drive operating configurations can be easily reproduced. Additionally, 
since none of these models contain inertia terms, simple algebraic equations are adequate to 
fully define model performance. Based on the results from my investigation of these five 

models, the applicability of different harmonic-drive modeling strategies will now be 
summarized. 

I began my investigation with an analysis of an ideal transmission model. Not 
surprisingly, since this representation failed to capture any of the frictional losses of the 
actual system, predicted operating velocities were orders-of-magnitude higher than in 
reality. To improve performance, I then included frictional losses into the dynamic 
representation. Specifically, dynamic friction, velocity-independent friction, and position- 
dependent friction were all appended to the ideal transmission model. As expected, this 
model predicted velocities which matched experimental observations in most ranges but 
failed to capture the resonant vibration observed in the actual system. Based on these 
results, it was concluded that, for systems where resonant behavior does not appear within 
the operating envelope of the harmonic drive, this simple representation may suffice. 

Unfortunately, in the three transmissions studied, system resonance appeared 
frequently and dominated behavior over a substantial range of harmonic-drive operating 
velocities. In order to capture these effects, both non-linear compliance and kinematic error 
were incorporated into the system model. By using stiffness values estimated from 
dynamic observations, the resulting predictions of resonant behavior from the model agreed 
closely with the experimental results. Unfortunately, since this improved harmonic-drive 
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model was incapable of describing the losses which occur at system resonance, model 
predictions at velocities which excited resonance vibration were still inaccurate. Based on 
these results, it was concluded that this model provides a slightly more accurate 
representation than it predecessor for describing the dynamic behavior of systems which 
experience little resonance vibration. 

In a final attempt to improve model performance, a representation was developed 
which incorporated a simplified gear-tooth meshing mechanism with Coulomb friction at 
the gear-tooth interface. Because of enhanced dissipation from the Coulomb friction at 
system resonance, this model better reproduced the actual velocities and torques that 
occurred in resonant regimes. In particular, the sudden velocity jumps that were observed 
in the experimental response when velocities exceeded resonant regions also appeared in 
the simulated response. From these optimistic results, it was concluded that the gear-tooth 
meshing mechanism plays a significant role in harmonic-drive operation and cannot be 
ignored in an accurate harmonic-drive model. Unfortunately, since some torque levels and 
velocity fluctuations predicted by this model differed from experimental observations, the 
need for more detailed representations was demonstrated. Nevertheless, it was concluded 
that, with careful estimation of input parameters, this model may prove useful in 
applications in which understanding harmonic-drive operation through resonant regions is 
required. 

From the results of the five transmission models, the complexity of the task of 
harmonic-drive modeling has been exposed. Not only are friction, compliance, and 
kinematic error essential factors for describing harmonic-drive behavior, but the gear-tooth 
meshing dynamics also play an essential role in determining the overall performance of the 
transmission. Additionally, experimental results indicate that the actual values of these 
transmission properties can vary significantly with environmental and operating conditions, 
which makes accurate characterization of these effects difficult if not impractical. To make 
matters worse, not only are reliable measurements of transmission properties difficult to 
obtain, but, for stiffness and coulomb friction in particular, precise estimates are critical for 
good model performance. In fact, the only way in which accurate values could be 
determined for the parameters in this model was by measuring carefully the dynamic 
response over a wide range of operation and adjusting of model parameters to match this 
observed response. This unfortunate result reveals that reliable harmonic drive models 
cannot be developed unless experimental data from the actual harmonic-drive system are 
available. In the face of these modeling obstacles, it seems unlikely that any model which 
does not capture the physics of harmonic-drive operation to a much greater level of detail 
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will ever be able to deliver accurate and reliable performance over all operating regions. 
Given this conclusion, I hope the models developed in this investigation will act as a 
springboard from which to launch more detailed future explorations. 


4.3 Recommendations 

Given the conclusions of my experimental and theoretical investigation, I can offer 
a few recommendations to designers wishing to model systems employing harmonic 
drives. The first step in understanding the extent to which harmonic-drive dynamics may 
influence system performance is it to determine the prominence of resonance vibration 
within the operating range of the transmission. Due to the unreliability of static 
measurements and catalog predictions of stiffness values, resonance vibration can be 
accurately located only during actual system operation. However, an order-of-magnitude 
estimate for the natural frequency of the system can be made by using linearized catalog 
stiffness values and system inertia measurements. If a harmonic of the rotational velocity 
of the wave-generator relative to the flexspline coincides with the estimated natural 
frequency of the system at any point in the operating envelope of the transmission, it is 
likely that resonant vibration will occur. Since resonant vibration greatly influences 
operating behavior and is difficult to model, the best solution to this dilemma is to revise 
the machine specifications to avoid these ranges of operation. 

For researchers embarking upon experimental investigations of harmonic-drive 
behavior, I can make a few suggestions about worthy areas for exploration. Specifically, 
in order to provide accurate stiffness data for harmonic-dnve models, a reliable method for 
removing the friction component from static stiffness measurements would be enormously 
valuable. Alternatively, modeling efforts would also benefit from an experimental 
procedure that could reliably measure harmonic-drive compliance dynamically. Since 
harmonic-drive stiffness and friction can vary considerably with such factors as gear-tooth 
preload and wave-generator angle, empirical data describing these variations could be very 
beneficial for anticipating ranges of parameter variation in different harmonic drives. 
Additionally, since many harmonic-drive properties are influenced by rotational direction, 
an experimental investigation in this area could also prove worthwhile. Lastly, since 
dynamic representations have been shown to depend strongly on accurate friction 
measurements, a method for experimentally separating different friction components 
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present in harmonic drives would be useful for improving performance of harmonic-drive 
models. 

Lastly, for future harmonic-drive modelers, I propose a few areas for investigation 
that seem most fruitful. In particular, as described in detail previously, a harmonic-drive 
representation that can describe the transition that occurs when the torque at the gear teeth 
changes direction should show improved results over the models in this document. After 
this improvement has been implemented, a model which predicts the influence of gear-tooth 
preload on transmission friction and compliance is also likely to improve dynamic 
response. Finally, since the kinematic error in harmonic-drives probably influences the 
gear-tooth preload itself, a representation which incorporates this behavior will more 
accurately represent the actual physical mechanism inside the harmonic-drive transmission. 
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This section contains citations for literature referenced in this thesis as well as 
citations for related harmonic-drive research. For convenience, the list of references on 
harmonic-drive research has been subdivided into categories describing the nature of the 
literature and popular areas of study. Papers which could fit reasonably into more than one 
category are listed only once. 

The references presented below were compiled from International CD-ROM 
database searches for harmonic-drive-related literature published in the past ten years. 
Additional references were obtained from library literature searches and reference lists from 
individual papers. The majority, but not all of these papers were located in the MIT 
Libraries. 
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Appendix A: Sensor Calibration 


This appendix describes the calibration procedures and results for all of the 
experimental sensors. 


A.1 Encoders 

Using processing boards to count encoder pulses, encoder measurements were 
directly translated into degrees of rotation. In order to make this simple conversion, the 
Whedco boards were initialized with the number of lines on each encoder wheel and the 
number of counts to extract for each line. These parameters are listed in table A. 1 for each 
encoder. The three input encoders directly measured motor rotation, while the output 
encoder on joint 1 monitored output rotation through an 8:1 gear reduction. 


Table A.1: Encoder < 

calibration parameters 



1 counts 



iiiiiiiiiiiii 

•:::->r*:-r;r>;r*r*i;r;r*:;.‘jr*r;r*.*^r*r*rfr*r*r*r*r*r*r*r*r*r*r*r*r*r»Sr*r*r*r*r*r*;f.**r*:*r*:'.**r*r*r*r*j;r;r*r*i*r»i*r*r*r*r*r’r*r-i*:-r*r-r* 



Joint 1 input encoder 

1000 

1,2, or 4 

«Joint i output encoder 

2500 

1,2, or 4 





1000 

1,2, or 4 

Joint 3 input encoder 

500 

1, 2, or 4 


A.2 Resolvers 

The resolver position information was returned as a 16-bit number representing a 
fraction of a complete rotation of the resolver shaft. Consequently, to convert this signal 
into degrees of resolver rotation, the number of resolver counts was multiplied by (360.0 
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degrees / 2 16 counts). The resulting value was then converted to degrees of output rotation 
by factoring in the sensor gear-ratio of 6.66, on joint 2, or 7.2, on joint 3. 

The resolver converter electronics could also return an analog voltage proportional 

to the resolver velocity. Since this information was never used, the resolver velocity was 
not calibrated. 


A.3 Tachometers 

The tachometers on joint 1 and 2 were calibrated by simultaneously collecting raw 
velocity data, measured in A to D counts, and encoder position data, measured in degrees. 
The tachometer data was then integrated to yield raw position data which was compared to 
the actual encoder position data to derive a calibration factor. Figures A.l and A.2 show 
the linear fits of integrated tachometer position data versus encoder data for the two 

tachometers. The resulting calibration factors and their associated uncertainty are 
summarized in table A.2. 
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Appendix A: Sensor Calibration 



Jolm 1 fftput tachometer 


Joint 2 inpat tachometer 


liilllil 

!!!!! 81 p!!lll 

16.85 

0.4 % 

17.33 

0.4 % 


A.4 Current Sensors 

The current produced by the amplifiers used to drive the DC motors was monitored 
by the A to D board. A to D counts were calibrated against amperes by recording the 
number of counts for a range of known currents measured with an ammeter. The resulting 
calibration values are summarized in table A.3 for the three joints. 
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Table A.3: Current sensor calibration factors 


IlIPiilliBI 



lllllilillllllll 






0.0115 





0.0122 

Joint 3 Current Censor 


0.0115 


A.5 Torque Sensors 

Torque sensor signals were processed by a conditioning circuit which returned a 
voltage proportional to the torque experienced by the sensor. This voltage was monitored 
by a 12-bit A to D channel which read the signal in A to D counts. This A to D reading was 
calibrated against actual torque measurements produced by a pre-calibrated, JR-3 six-axis 
force gauge. The calibration procedure for all three joints began by rigidly locking the 
wave-generator to the circular spline using a specially designed fixture. As illustrated in 
figure A.3, the force sensor was mounted on the end of the output link for each joint and a 
force was applied in the direction of joint rotation. By multiplying the force reading by the 
radius of the applied load, an accurate measurement of the resulting torque was found. 
Using this method, torque data was collected from both the JR-3 force sensor and the 
uncalibrated torque sensor over a range of applied loads in both directions. This data is 
shown in figures A.4 through A.9 for the three sensors in both loading directions. The 
calibration factors derived from averaging the linear-fits to this data in both direction are 
summarized in table A.4 for each testing station. 


I Table A.4: 

Torque sensor calibration factors 



H+m 

count 

IMS 

uncertainty 


count 

Joint 1 torque sensor 

0.0437 

0.3864 

0.8 % 

Joint 2 torque sensor 

0.0290 

0.2567 

0.8 % 

Joint 3 torque sensor 

0.0065 

0.0572 

1.8 % 

























































































































































Force 


Figure A.3: Experimental apparatus used to calibrate torque sensors 
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Figure A.7: Joint 2 torque sensor calibration data for negative loading 



Figure A.8: Joint 3 torque sensor calibration data for positive loading 
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Force Sensor Torque (in*lbs) 
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Appendix B: Simulation Code 


The computer code that was used to numerically solve the differential equations 
describing the dynamic behavior of the different harmonic-drive models is listed in this 
appendix. The header-file of variable definitions used by this program is also presented 
after the main program listing. In addition, this program relies on a function, not listed 
here, which can integrate ordinary differential equations. The function I used was a 
Runga-Kutta solver with adaptive stepsize control taken largely from reference ( 85 ), 
Numerical Recipes in C. During execution, this program reads the input parameters from a 
file, calculates and integrates the specified equations of motion, and outputs time-response 
data for the variable names specified in another file. I hope that this code can serve as 
template for rapid development of future harmonic-drive simulations. 


/ This is file hd_model.c. It contains the simulation code used to solve a 

* system of differential equations which characterize the behavior of a 
harmonic drive. Included in this file are two function which calculate 
equations of motion for the two different joint configurations used for 
the three harmonic-drive testing stations. The first two joints (shoudler 

* and elbow) use configurationl and the third joint (wrist) uses the second 
configuration. Additionally, code to calculate the equations for two 

fferent harmonic-drive models is also included in this file. This code 
is located in three separate functions which can be called arbitrarily 
by any function which calculated the equations of motion of a dynamic 

* system. If desired, this program can also return the energy distribution 
in the simulation and the stiffness of the given model. */ 

#include <stdio.h> 

♦include <stdlib.h> 

♦include <math.h> 

♦include <string.h> 

♦include "hd_model.h" 

/* Declare and initialize the global variables, */ 

int debug = FALSE; 

int calc_energy = FALSE; 

int print_energy = FALSE; 

int stiffness_flag = FALSE; 

int stiffness_flag2 = FALSE; 

int ideal_flag = FALSE; 

int fast_plot; 

int plot = FALSE; 

int model; 

int joint = JOINT TO TEST;_ 
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I int num_selected_output var = 0 ; 
int system_order = ORDER; 
int index[NUM_OF_OUTPUT VAR]; 
float current_time; 

char graph_title[NUM_OF_OUTPUT VAR][50]; 
char graph_units[NUM_OF_OUTPUT~VAR][50]; 
char output_datafile[NUM_OF_OUTPUT VAR][100]; 
char plot_datafile{NUM_OF_OUTPUT VAR] [100] ; 

/* Initialize some useful character strings. ★/ 

char gnuplot_filename[] = "plot_commands.gnu"; 

char *joint_name[NUM_OF_AXIS] = {"shld", "elb", "wrist"}; 

char out-datapath[] = "/home/me/tut/robot/hd/model/final_model/simulation data/"- 
char indatapath[] = "/home/me/tut/robot/hd/model/final_model/simulation code/"-' 
char plotdatapath[] = "/home/me/tut/robot/hd/model/final_model/plot_data/"; 

/* Initialize the array of maximum velocity contraints on the non-linear 
damping functions for each model. */ 

float max_f riction_vel [ HD_M0DELS ] [NLJM_OF_AXIS] = {{ 135.0, 300.0, 500.0}, 

(250.0, 550.0, 1000.0}}; 

'* Declare the functions which contain the equations of motion and the 
Runge-Kutta solver function which is defined in the file rkqc.c. */ 
foid derivs_configl(double, double *, double *); 
roid derivs_config2(double, double *, double *); 
ixtern void rkqc(double *,double *,int,double *,double, 

double,double *,double *,double *,void (*)()); 

Declare the functions which calculate the energy for both of the two 
* joint configurations and the total system energy. */ 

oid calculate_config_energy(double); 

oid calculate__total_energy (double *, double *); 

* Declare functions which calculate the total equations of motion for the 

* har | nonic <J rive system when the transmission is ideal. in this case, the 
system reduced to a first-order system. */ 
oid derivs_ideal_configl(double, double *, double *); 
oid derivs_ideal_config2(double, double *, double *); 

* * P ° inte \ to . a function that will hold the appropriate pointer to 

the function which will calculate the equations of motion. */ 
oid (*ptr_derivs)(double, double *, double *); 

* p D !n!t re the j- ndividual functions which calculate the input parameters, dynamic 
equations and energy values for the different harmonic-drive models. */ 

3id calculate_rotary_params(void) ; 

3id rotary_hd_model(double, double, double, double, double, double, 

double *, double *, double *); 

aid rotary_hd_energy(double); 

)id calculate_gear_tooth_params(void) ; 

.id gear_tooth_hd_model(double, double, double, double, double, double, 

double *, double *, double *); 

)id gear_tooth_hd_energy(double) ; 

N ° W lnit ialize these functions into arrays of function pointers for 
easy access. */ 

'id (*calculate_params(HD_MODELS])(void) 

•w 7*^ alculate - rotary _ pararns . calculate_gear tooth params}; 
id (*hd_model(HD_MODELS]l(double, double, double, double, double, double, 

double *, double *, double *) 

= {rotary_hd_model, gear_tooth_hd_model}; 
id (*hd_energy[HD_MODELS])(double) 

= {rotary_hd_energy, gear_tooth_hd_energy}; 
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/* Declare the functions which (1) calculate stiffness data, (2) model the 

* saturation limits of the current amplifiers, and (3) determine the static 

* and dynamic friction values, respectively. */ 
void collect_stiffness_data(void) ; 

double iregulate(double, double, double); 

double calculate_friction(double, double, double, double, double, 

double, double, double, double, double); 

/* Declare the rest of the functions. */ 

void get_arguments(int argc, char **argv); 

void usage(void); 

void read_input_data(FILE *) ; 

void read_input_filenames(FILE *); 

void print_input_parameters(void) ; 

void print_variables(void) ; 

void export_vectors(FILE *, double *, double *, int, char *); 
void generate_gnuplot_file(char *) ; 


/* 

, */ 

/* Main Program 

/* 

, V 

/★*****************************», 

/* This main program opens and reads the input data files, initializes appropriate 
* variables and prints/plots the results. */ 
main(int argc, char **argv) 

{ 

int h, i, j; 

int num_of_trials[NUM_OF_AXIS] = {SHLD_NUM_OF_DATA_RUNS, 

ELB_NUM_OF_DATA_RUNS, 

WRIST_NUM_OF_DATA_RUNS}; 
int good_steps[1], bad_steps[1] ; 

float input_amps [ N(JM_OF_AXIS ] [ SHLD_N(JM_OF_DATA_RUNS ] = 

{{1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.4, 5.0, 6.0}, 

{1.4, 1.8, 2.2, 2.6, 3.0, 3.4, 3.8, 4.5, 0.0, 0.0}, 

{0.28, 0.36, 0.40, 0.44, 0.48, 0.52, 0.56, 0.60, 0.0, 0.0}}; 
char unix_command[200] ; 
char vector_name[100] ; 

char *input_datafile[NUM_OF_AXIS J = {"hd_shld_input_data.dat", 

"hd_elb_input_data.dat", 

"hd_wrist_input data.dat"}; 

FILE *infile; 

FILE *outfile[NUM_OF_OUTPUT_VAR]; 
extern int kount; 


/* Read and process the command-line arguments. */ 
get_arguments(argc, argv); 

/* Open the input data file for the given joint, read the simulation 
parameters from the file, and then close the input file. */ 
open_file(& in file, indatapath, input_datafile[joint]); 
read_input_data(infile); 
close_file(infile); 

/ Calculate necessary parameters for the given model from the values read 
* from the input file. */ 

(*(calculate_params[model)))(); 

/* Print out the input data if in debugging mode. */ 
if (debug) 
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print_input_parameters {) ; " " ~ —-- 

/* Set the pointer to the appropriate function containing the equations of 
* motion. */ 

if ( joint == WRIST) 

{ 

if(ideal_flag) 

ptr_derivs = derivs_ideal_config2; 
else 

ptr_derivs = derivs_config2; 


else 


if (ideal_flag) 

ptr_derivs = derivs_ideal_configl; 
else 

Ptr_derivs = derivs_configl; 


{ If stiffness data is to be collected, jump to the appropriate function. */ 
if (stiffness_flag) 

collect_stif fness_data(); 

/* N ° W rSad ln the data from the other input file that provides information 
about what data to output. */ 

° P en_f ile (4 in file, indatapath, "hd_input_f ilenames .daf) ; 
read_input_filenames(infile); 
close_file{infile); 

/* Now if the largest line number in the index(] vector is greater than 
* ° r equal to the first line of the energy variables in the second 
input data file, then energy is required as an output, so that the 
energy flag should be set so that energy is calculated. */ 
if (index [ num__selected_output_var-1 ] >= TOTAL ENERGY) 
calc_energy = TRUE; 

7 * ^ 1S tC> bS calculat ed, increase the order of the system 

to ORDER+1 in order to accomodate the additional total energy variable 
which is to be integrated as a member of the state variable. */ 
if(calc_energy) 

system_order = ORDER + 1 ; 

/* If the -plot option is not selected, set the number of trials to 1 */ 

if ( !plot) ’ 

num_of_trials[joint] = 1 ; 

/* N ° w . run the simulation num_of_trials-times by calling odeint(). Function 
* ?' ^° ntains the R unga-Kutta numerical solver which calls the functions, 

listed below which calculate the equations of motion for the selected dynamic 

f °r(j = 0 ; j < num_of_trials[joint]; j++) 

/* Change the value of the requested motor current if necessary */ 
if (plot) ' 

irequested = input_amps[joint][j]; 

/* Restore or initialize the xstart[) values. */ 
for(i = 0; i < (ORDER+1); i++) 
xstart[i] = xstart_sav[i]; 

printf("\n Starting simulation at %f amps- \n", irequested); 

odeint(xstart, system_order, initial_time, final_time, EPS, INITIAL STEP 
MINIMUM_STEP, good_steps, bad_steps, ptr_derivs, rkqc); 
printf("\n Simulation complete. \n"); 
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if(debug) 

{ - 

printf("\n Number of good steps taken: %d", good_steps[0]); 
printf ("\n Number of bad steps taken: %d\n", bad steps[0)); 


/* Initialize all output files, export the desired output data to them, 

* and close the files. If the -plot option is selected, then save 
the specified data in the plot_data directory with the appropriate 

* filename. */ 

for(i =0/ i < num_selected_output var; i++) 

{ 

if(plot) 


sprintf(plot_datafile ( i], "%s.%4.2famps", 

output_datafile[i), input_amps[joint][j])/ 
create_file(&(outfile[i]), plotdatapath, plot_dataf 
export_vectors(outfile[i], time, output_data[i+1], 
output_datafile[i]) ; 


ile[i[ 
kount, 


else 

{ 

create_file(&(outfile[i]), outdatapath, outpu 
export_vectors(outfile[i], time, output_data[ 
output_datafile[i]) ; 

) 

close_file(out file[i]); 


t_data 
i + 1 ] , 


file [. 
kount. 


/* Now plot the results using xgraph, if slow printable plotting is desired, 
* or gnuplot, if the fast_plot option is selected. */ 
if((plot) 

{ 

if(fast_plot) 

{ 

generate_gnuplot_file(gnuplot_filename); 

sprintf(unix_command, "gnuplot %s%s", outdatapath, gnuplot_filename); 
printf <"\n Executing the Unix command: %s\n", unix command); 
system(unix_command); 

) 

else 

{ 

for(i = 0; i < num_selected_output_var; i++) 

{ 

sprintf(unix_command, "xgraph %s%s -t Harmonic_Drive_Simulation\ 
-y %s_%s -x Time_sec &", outdatapath, output_datafile[i], 
graph_title[i], graph_units[i] ) ; 
printf ("\n Executing the unix command: %s\n", unix command); 
system(unix command); 




**★★★***■★*•**★*•*■ 


/* Functions which describe the dynamics of the two joint configurations, 
/* 


*■*•*•*■*★★★*■*****•***•*** 
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/* Function derivs_configl{) can be used to solve the equations of motion for the 
harmonic drive joint configuration in which: (1) the motor drives the wave- 
generator and is mounted to the circular spline, (2) the flexspline is attached 
to ground, and (3) the circular spline is the output port. */ 
void derivs_configl(t, x, dxdt) 
double t; 

double x[], dxdt[]; 

{ 

/* Set the global time variable to the current time. */ 
current_time = t; 

/* Set the local variable definitions to the value of the state variable 
* determined on the previous step. */ 

v[POS_IN] = x[POS_IN_INDEX] ; 

v [ POS_OtJT ] = x [ POS_OUT_INDEX ] ; 

v[VEL_IN] = x[VEL_IN_INDEX] ; 
v[VEL_OUT] = x[VEL OUT INDEX]; 


/* From these positions and velocities, call the appropriate harmonic drive 
* model to calculate the resulting torques. */ 

(*(hd_model[model]))(v[POS_IN], 0.0, v[POS OUT], 

v[VEL_IN], 0.0, v[VEL_OUT], 

&(v[TORQUE_WG]), &{v[TORQUE_FS]), &(v[TORQUE CS]) ) ; 


/* Set the value of torque sensor to the torque read 
v[TORQUE_SENSOR] = -v[TORQUE FS]; 


on the flexspline. 


'/ 


/* Calculate the motor torque as a function of current: 

Note that the current, factual, is calculated by the function 
iregulate() which takes into account the non-linear behavior of the 
current amps. Note that v[CURRENT_SENSOR] must be initialized to zero 
Recall that the relative velocity of the motor is the absolute motor 
velocity minus the absolute output velocity. */ 
v[CURRENT_SENSOR] = iregulate(v[CURRENT_SENSOR], irequested, 

(v[VEL_IN] - v[VEL_OUT])); 
v[TORQUE_MOTOR] = motor_kt * v[CURRENT_SENSOR]; 

/ Calculate the input and output damping. */ 
v[TORQUE_B_IN] = b_in * (v[VEL_IN] - v[VEL OUT]); 
v[TORQUE_B_OUT] = b out * v[VEL OUT]• 


/ Finally, calculate the equations of motion. */ 
dxdt[POS_IN_INDEX] = v[VEL_IN]; 

dxdt[VEL_IN_INDEX] = (1.0 / (CONV * inertia_in)) * (v[TORQUE MOTOR] - 

v[TORQUE_B_IN] - 

, v[TORQUE WG]); 

dxdt[POS_OUT_INDEX] = v[VEL_OUT]; 

dxdt[VEL_OUT_INDEX] = (1.0 / (CONV * inertia_out)) * (v(TORQUE CSJ - 

v[ torque_b_out] + 
v[TORQUE_B_IN]); 

/ Now calculate the values that the sensors actually read. */ 
v[INPUT_POSITION_SENSOR] = v[POS IN] - v[POS OUT]; 
v[INPUT_VELOCITY_SENSOR] = v[VEL IN] - v[VEL OUT]; 
v[OUTPUT_POSITION_SENSOR] =v[POS_OUT); 

V[OUTPUT_VELOCITY_SENSOR] = v[VEL~OUT]; 

v[POSITION_ERROR] = (V[INPUT_POSITION_SENSOR] / N) - v[OUTPUT_POSITION SENSOR] 

/* If the total energy is selected as an output variable, call 
the function that calculates the total power and integrates 
* it using the runga-kutta solver. */ 
if(calc_energy) 

calculate tota l energy(x, dxdt); 
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/* If the debugging option is selected, call the function to print 
* out some variables. */ 
if(debug) 

print_variables(); 


/* Function derivs_config2() can be used to solve the equations of motion for 

* the harmonic drive joint configuration in which: (1) the motor drives the 

* wave-generator and is mounted to the circular spline, (2) the flexspline is 

* attached to the output link, and (3) the circular spline attached to ground. */ 
void derivs_config2(t, x, dxdt) 

double t; 

double x[], dxdt[]; 

{ 

/* Set the global time variable to the current time. */ 
current_time = t; 

/* Set the local variable definitions to the value of the state variable 

* determined on the previous step. */ 

V[POS_IN] = x[POS_IN_INDEX]; 

v[POS_OUT] = x(POS_OUT_INDEX] ; 

v[VEL_IN] = x[VEL_IN_INDEX] ; 

v[VEL_OUT] = x[VEL_OUT_INDEX] ; 

/* From these positions and velocities, call the appropriate harmonic drive 

* model to calculate the resulting torques. */ 

(*(hd_model[model]))(v(POS_IN], v[POS_OUT], 0.0, 

v[VEL_IN], v[VEL_OUT], 0.0, 

&(v[TORQUE_WG)) , &(v{TORQUE_FS]), &(v[TORQUE_CS)) ) ; 

/* Calculate the motor torque as a function of current: 

* Note that the current, iactual, is calculated by the function 

* iregulateO which takes into account the non-linear behavior of the 

* current amps. Note that v[CURRENT_SENSOR] must be initialized to zero. */ 
v[CURRENT_SENSOR] = iregulate(v[CURRENT_SENSOR], irequested, v[VEL IN]); 

v[TORQUE_MOTOR] = motor_kt * v(CURRENT_SENSOR]; 

/* Calculate the input and output damping. */ 
v[TORQUE_B_IN] = b_in * v(VEL_IN]; 
v[TORQUE_B_OUT] = b_out * v[VEL_OUT); 

/* The torque sensor is mounted between the motor and circular spline and 

* ground, and therefore sees the reaction torque from the motor, the motor 

* damping torque, and the circular spline torque. */ 

v [TORQUE_SENSOR] = v[TORQUE_CSj - v ( TORQUE_MOTOR ] + v [ TORQUE^_B_IN ] ; 

/* Finally, calculate the equations of motion. */ 
dxdt[POS_IN_INDEX] = v[VEL_IN]; 

dxdt[VEL_IN_INDEX] = (1.0 / (CONV * inertia_in)) * (v[TORQUE MOTOR] - 

v[torque_b_in7 - 

v[TORQ(JE_WG]); 

dxdt[POS_OUT_INDEX] = v[VEL_OUT]; 

dxdt [ VEL_0(JT_INDEX ] = (1.0 / (CONV * inertia_out) ) * (-v[TORQUE FS] - 

v[TORQUE_B_OUT]); 

/* Now calculate the values that the sensors actually read. */ 
v[INPUT_POSITION_SENSOR] = v(POS_IN]; 
v[INPUT_VELOCITY_SENSOR] = v[VEL_IN]; 
v[OUTPUT_POSITION_SENSOR] = v[POS_OUT]; 
v[OUTPUT_VELOCITY_SENSOR] = v[VEL_OUT]; 

- ^XPQS^TI0N^^RR0RJ__=^_^vJ^NPUT_P 0S_IT_I0N__ SENSOR ] / N) -t- v [OUTPUT POSITION SENSOR]; 
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/* If the total energy is selected as an output variable, call 
the function that calculates the total power and integrates 

* it using the runga-kutta solver. */ 
if(calc_energy) 

calculate_total_energy(x, dxdt); 

/* If the debugging option is selected, call the function to print 

* out some variables. */ 
if(debug) 

print_variables (); 


*********** ******** 


/****************************_^^^^^^^ 

/* 

/* Functions which describe the energy of the two joint configurati 
/*******************.**.*.***.**** #A ^^ 


********************* ********/ 

*/ 

ions. */ 


*****************, 


****************************^*^^ 


/* F ^ nctlon , calcul ate_total_energy<) calculates the total power in the system 

* ° f ^ uatlons Presented in derivs_configl() or derivs_config2(). This power 

* de ^rmined on every time step and stored in the state vector which is 

* ^° 1Ved by the Run 9e-Kutta solver. This can be used to verify that energy 
is being conserved in both dynamic models. */ 

void calculate_total_energy(x, dxdt) 
double x[], dxdt(]; 

{ 

double power_b_in, power_hd_out; 

/* First store the integrated power in the total energy variable */ 
v [ TOTAL_ENERGY ] = x [ENERGY]/ vanaoie. / 

/* Now calculate the input damping power loss depending on the joint 

* conFl ^ration. Also calculate the appropriate torque on the output 
port of the harmonic drive. */ 
if ( joint == WRIST) 

{ 

power_b_in = v[TORQUE_B_IN] * (v[VEL_IN]*CONV)/ 
power_hd_out = -v[TORQUE_FS] * (v[VEL OUT]*CONV); 


power_b_in - v[TORQUE_B_IN] * ((v[VEL_IN] - v[VEL OUT])*CONV)• 

^ power_hd_out = v[TORQUE_CS] * (v(VEL_OUT]*CONV) ; _ 

/* Calculate the total instantaneous power in the system. */ 
dxdt[ENERGY] = ((v[TORQUE_MOTORJ * (v[VEL IN]*CONV)) - 

(inertia_in * (v[VEL_IN]*CONV) * (dxdt[VEL_IN_INDEX] *CONV)) - 
(inertia_out * (v[VEL_OUT]*CONV) * (dxdt[VELOUT INDEX ]*CONV)) - 
(v[TORQUE_WG] * (v[VEL_IN]*CONV)) + 

(power_hd_out) - 
(power_b_in) - 

(V[TORQUE_B_OUT] * (v[VEL_OUTJ*CONV))); 


/ l I^h ti0n ^ alculate - confi 9- ene rgy() is called by the function odeintO after 
* ^ PlSte tlme Ste ^ This function takes the current simulation time 

step_an d uses a simp le Euler integration scheme to determine the 


enerc 


in 
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* the system, excluding the harmonic drive. Since this function relies on 

* values calculated in the harmonic drive energy function, it should be called 

* after the specific harmonic-drive energy function is called. Due to the 

* similarities between the equations in derivs_configl() and derivs config2(), 

* this energy function is valid for both models. */ 
void calculate_config_energy(delta_t) 

double delta_t; 

{ 

double vel_in_rad, vel_out_rad; 

/* Be sure that all velocities and positions are in radians to ensure 

* proper unit-matching. */ 
vel_in_rad = v[VEL_IN] * CONV; 
vel_out_rad = v[VEL_OUT] * CONV; 

/* Calculate the energy input by the motor to the system by integrating 

* its power input over time. */ 

v[MOTOR_ENERGY] = v[MOTOR_ENERGY] + (delta_t * v[TORQUE_MOTOR] * vel_in rad); 

/* Calculate the kinetic energy of the inertias from (I *(w~2)/2) . */ 
v[KINETIC_ENERGY_IN] = (inertia_in * vel_in_rad * vel_in_rad) / 2.0; 
v [ KINETIC_ENERGY_OCJT ] = (inertia_out * vel_out_rad * vel_out_rad) / 2.0; 

/* Calculate the damping energy by integrating the power dissipation over 

* time. */ 

if(joint == WRIST) 

v[INPUT_DAMPING_ENERGY] = (v[INPUT_DAMPING_ENERGY] + 

(delta_t * (v[TORQUE_B_IN] * vel_in rad))); 

else 

v [ INPUT_DAMPING_ENERGY ] = (v [ INPUT__DAMP ING_ENERGY] + 

(delta_t * (v[TORQUE_B_IN] * 

(vel_in_rad - vel_out_rad)))); 
v[OUTPUT_DAMPING_ENERGY] = (v[OUTPUT_DAMPING_ENERGY] + 

(delta_t * (v[TORQUE_B_OUT] * vel_out_rad))); 

/* Now calculate some aggregate quantities. */ 

v[TOTAL_KINETIC_ENERGY] = v[KINETIC_ENERGY_IN] + v(KINETIC_ENERGY_OUT]; 
v[TOTAL_INPUT_ENERGY] = v[MOTOR_ENERGY]; 
v[TOTAL_LOST_ENERGY] = (v[INPUT_DAMPING_ENERGY] + 

v[WG_FRICTION_ENERGY j + 
v[TOOTH_TIP_FRICTION_ENERGY] + 
v[TOOTH_SURFACE_TOTAL_LOSS_ENERGY] + 
v[OUTPUT_DAMPING_ENERGY]); 

/* Now print out the energy values if the print_energy option is selected. */ 
if (print_energy) 

{ 

printf("\n total energy: %25.201f", v[TOTAL_ENERGY]) ; 

printf("\n total kinetic energy: %20.151f", v[TOTAL_KINETIC_ENERGY]); 

print f("\n total potential energy: %20.151f", v[TOTAL_POTENTIAL_ENERGY]) ; 
printf("\n total input energy: %20.151f", v[TOTAL_INPUT_ENERGYj); 

printf("\n total lost energy: %20,151f", v[TOTAL_LOST_ENERGY]); 

printf("\n motor energy: %20.151f", v[MOTOR_ENERGY]); 

printf("\n kinetic energy in: %20.151f", v[KINETIC_ENERGY_INj); 

printf("\n kinetic energy out: %20.151f'\ v[KINETIC_ENERGY_OUT]); 

printf("\n spring energy: %20.151f", v[SPRING_ENERGY)); 

printf("\n cyclic torque energy: %20.151f", v[CYCLIC_ENERGY]); 

printf("\n delta_t: %20.151f", delta_t); 

printf("\n"); 

} 

) 
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/* 


************* ******************* ir< . 1tlrif 


*********** 


/* Functions which describe behavior of system with ideal tr 
/* 

/*************************^.^___^_____^^^^^ 


ansmiss ion 


•AT*************. 


f * nnnf tl0n f eriv s_ldeal_config1() calculates the equations of motion for a 
configuration-1 sytstem that has an ideal harmonic-drive transmission. ■/ 
void derivs_ideal_configl(t, x, dxdt) 
double t; 

double x[], dxdt[]; 

{ 

/* Set the global time variable to the current time */ 
current_time = t; 

/ * Set the local variable definitions to the value of the state variable 
determined on the previous step. */ 

v[POS_IN] = x[POS_IN_INDEX]; 

v[POS_OUT] = x[POS_OUT_INDEX]; 

v[VEL_IN] = x[VE L_IN_INDEX); 

v[VEL_OUT] = x[VEL_OUT_INDEX]; 

/* Calculate the motor torque as a function of current. For the ideal 
model, assume that the current amps function ideally as well. */ 
v[CURRENT_SENSOR] = irequested; 
v[TORQUE_MOTOR] = motor_kt * v[CURRENT_SENSOR]; 

/* Calculate the input and output damping. */ 
v [TORQUE__B_IN ] = b_in * (v[VEL_IN) - vfVEL OUT]); 
v[TORQUE_B_OUT] = b_out * v[VEL OUT]; 

/* Finally, calculate the equations of motion */ 
dxdt[POS_IN_INDEX] = v[VEL_IN]; 

dxdt[VEL_IN_ IN DEX] = ((1.0 / (CONV * (inertia_in + 

(inertia_out/((N+1.0)*(N+1.0)))))) » 

(v[TORQUE_MOTOR] - v[TORQUE_B IN] - 

((v[TORQUE_B_OUT] - v[T0RQUE B IN])/(N+1.0) ) ) ) • 
dxdt[POS_OUT_INDEX] = dxdt[POS_IN_INDEX]/(N+l.0); 
dxdt[VEL_OUT_INDEX] = dxdt[VEL_IN_INDEX]/(N+l.0); 

/ Now calculate the values that the sensors actually read. */ 
v[INPUT_POSITION_SENSOR] = v[POS IN] - v( POS OUT]; 
v[INPUT_VELOCITY_SENSOR] = v[VEL IN] - v[VEL OUT]- 
v[OUTPUT_POSITION_SENSOR] = v[POS OUT]; 
v[OUTPUT_VELOCITY_SENSOR] = v[VEL OUT]; 

viTORQUE^S^- 0 ^^/ ^ VtI ^ UT - POSI ^ ION - SENSOR) / N) - v(OUTPUT POSITION SENSOR]; 
v IORQUE_FS] - (-N/(N + 1.0)) * (((CONV * inertia_in) * 

(dxdt[VEL_IN_INDEX]/(N+l.0))) + 

v [T OR QUE_B_OUT] - v[TORQUE B IN]); 
v[TORQUE_SENSOR] = -v[TORQUE FS]; — — 

/ trp^i^ Ulate thS torc ^ ues on the remaining harmonic-drive ports */ 
v[TORQUE_WG] = (-1.0/N) * v[TORQUE_FS]; P 

v [ TORQUE__CS ] = (-(N + 1.0) / N) * v[TORQUE FS ] ; 

/ \ ?! t0tal ener< ?Y is selected as an output variable, call 

the function that calculates the total power and integrates 
it using the runga—kutta solver. */ 
if <calc_energy) 

calculate_total_energy(x, dxdt); 

/* If thS debugging option is selected, call the function to print 
out some variables. */ 

if(debug) 
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print_variables() ; 

} 


/* Function derivs_ideal_config2() calculates the equations of motion for a 
* configuration-2 sytstem that has an ideal harmonic-drive transmission. */ 
void derivs_ideal_config2(t, x, dxdt) 
double t; 

double x[], dxdt[]; 

{ 

/* Set the global time variable to the current time. */ 
current_time = t; 

/* Set the local variable definitions to the value of the state variable 

* determined on the previous step. */ 

V[POS_IN] = x[POS_IN_INDEX); 

v[POS_OUT] = x[POS_OUT_INDEX]; 
v(VEL_IN] = x[VEL_IN_INDEX] ; 

v[VEL_OUT] = x[VEL_OUT_INDEX]; 

/* Calculate the motor torque as a function of current. For the ideal 

* model, assume that the current amps function ideally as well. * / 
v[CURRENT_SENSOR] = irequested; 

v[TORQUE_MOTOR] = motor_kt * v[CURRENT_SENSOR]; 

/* Calculate the input and output damping. */ 
v[TORQUE_B_IN] = b_in * v[VEL_IN]; 
v[TORQUE_B_OUT] = b_out * v[VEL_OUT]/ 

/* Finally, calculate the equations of motion. */ 
dxdt[POS_IN_INDEX] = v[VEL_IN]; 

dxdt[VEL_IN_INDEX] = ((1.0 / (CONV * (inertia_in - (inertia out/(N*N))))) * 

(v[TORQUE_MOTORJ - v[TORQUE_B_IN] - 
(v[TORQUE_B_OUT]/N))); 

dxdt[POS_OUT_INDEX] = dxdt[POS_IN_INDEX] / (-N); 
dxdt[VEL_OUT_INDEX] = dxdt[VEL_IN_INDEX] / (-N); 

/* Now calculate the values that the sensors actually read. */ 
v[INPUT_POSITION_SENSOR] = v(POS_IN] - v[POS_OUT]; 
v[INPUT_VELOCITY_SENSOR] = v(VEL_IN] - v[VEL~0UT]; 
v[OUTPUT_POSITION_SENSORJ = v[P0S_0UT]; 
v[OUTPUT_VELOCITY_SENSOR] = v[VEL_0UT); 

V[POSITION_ERROR] = (v[INPUT_POSITION_SENSOR] / N) + v[OUTPUT_POSITION SENSOR]; 
v[T0RQUE_FS] = -(((CONV * inertia_out) * (dxdt[VEL_IN_INDEX]/(-N))) t 
v[TORQUE_B_OUT]); 

/* Calculate the torques on the remaining harmonic-drive ports. */ 
v [TORQ(JE_WG ] = (-1.0/N) * v [ TORQUE_FS ] ; 
v[TORQUE_CS] = (-(N+l.0)/N) * v[T0RQUE_FS]; 

/* The torque sensor is mounted between the motor and circular spline and 

* ground, and therefore sees the reaction torque from the motor, the motor 

* damping torque, and the circular spline torque. */ 

V[TORQUE_SENSOR] = v[TORQUE_CS] - v[TORQUE_MOTOR] + v[TORQUE_B_IN]; 

/* If the total energy is selected as an output variable, call 

* the function that calculates the total power and integrates 

* it using the runga-kutta solver. */ 
if (calc_energy) 

calculate_total_energy(x, dxdt); 

/* If the debugging option is selected, call the function to print 

* out some variables. */ 
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if(debug) 

print_variables(); 


/************************************************************ 
/* 

/* Functions determining equations for two different 
/ * 

/*****************************************************, 


******************i 


harmonic-drive models 


******************** 


****************, 


/* Function rotary_hd_model() calculates the three-port harmonic-drive 
equations for a non-ideal rotational model. */ 
void rotary_hd_model(pos_wg, pos_fs, pos cs, 

vel_wg, vel_fs, vel_cs, 

torque_wg, torque_fs, torque cs) 
double pos_wg, pos_fs, pos_cs; 
double vel_wg, vel_fs, vel cs; 
double *torque_wg, *torque~fs, *torque cs- 
{ 

/* Store the velocities and positions and the fs, cs, and wg ports 
in the global data vector to begin calculations. */ 
v[POS_WG] = pos_wg; 
v[POS_FS] = pos_fs; 
v[POS_CS] = pos cs; 
v [ VE L_WG ] = ve 1 wg ; 
v[VEL_FS] = vel~fs; 
v[VEL_CS] = vel_cs; 

/* Calculate the continuity constraint of the cyclic torque */ 
v[POS_K] = v[POS_WG]; 4 

v[VEL_K] = v[VEL_WG]; 

/* Calculate some other continuity constraints. */ 
v(POS_N_FS] = v[POS_FS]; 
v[VEL_N_FS] = v(VEL_FS]/ 
v[POS_N_CS] = v[POS_CS]; 
v[VEL_N_CS] = v[VEL_CS]; 

thS kinem « ic funct ion and its derivative such that 

(d/dt)v[ERFN] = v[VEL_WG] v[DERF]. */ 

v[ERFN] = ( (error_amplitudeO * Sin(((1.0 * v[POS_WG]) + error phaseO) * CONV) ) • 

(error_amplitudel * Sin(((2.0 * v[P0S_WGj) + error_phasel) * CONV)) + 

vfDERFl - ( ( n r n° r *- a S u it * Ude2 * si n(((4.0 * v[POS_WG]) + er ror_phase2) * CONV) ) ) ; 
v [DERF ] - (d.o * CONV * error_amplitudeO * 

Cos(((1.0 * v[POS_WG]) + error_phaseO) * CONV)) + 

(2.0 * CONV * error_amplitudel * 

Cos (((2.0 * v(POS_WG]) + error_phasel) * CONV)) + 

(4.0 * CONV * error_amplitude2 * 

Cos (((4.0 * v[POS_WGJ) + error_phase2) * CONV))); 

/* Now calculate the position and velocity equations imposed by the 

vrPOsTwGi r - M7!T 1C 7 d m Ve gear " reduction with Position error included. */ 
v[POS_N_WG] - (((N + 1.0) * v[POS N CSJ) - 

(N * v[POS_N_FS]) + 

(N * V(ERFN] ) ) ; 

v[VEL_N_WG) = (((N + 1.0) * v[VEL N CS]) - 
(N * v[VEL_N_FS]) + 

(N * v[VEL_WG] * v[DERF])); 

^ aiCudate ^atic, dynamic, and cyclic friction which acts between the 

_ circular spline and flexspline */ 
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v[VEL_KD_FRICTION] = v[VEL_N_CS] - v[VEL_N_FS]; 

V[TORQUE_HD_FRICTION} = 

calculate_friction(v[VEL_HD_FRICTION], 

b_hd_constant, b_hdl, b_hd2, 
stiction_torque_hd, stiction_vel_hd, 
v[POS_OUT], 

cyclic_f riction_amp_hd, cyclic_f riction_phase_nd, 

(max_friction_vel[model][joint])); 

/* For the sake of continuity in the simulation, ramp the constant friction 

* in the harmonic drive from zero to the desired value if at the 

* beginning of the step-response trial. */ 
if ( ! (stif fness_f lag || stif f nessj lag2) ) 

if (current_time < 0.05) 

v[TORQUE_HD_FRICTION] = ((((current_time - 0.05)/0.05) + 1.0) * 

v[TORQUE_HD_FRICTION]); 

/* Find the harmonic drive stiffness coefficients: */ 
v[Kl] = (kl_constant + 

kl_amplitude * Sin(((2.0 * v[POS_WG]) + kl_phase) * CONV)); 
v[K2] = (k2_constant + 

k2_amplitude * Sin(((2.0 * v(POS_WG)) + k2_phase) * CONV)); 

/* Calculate the torque on the spring. */ 
v[TORQUE_K] = ((v[K1] * (v[POS_K] - v[POS_N_WG])) + 

(v[K2] * Power((v[POS_K] - v[POS_N_WG]), 3)))/ 

/* Calculate the compatibility, or torque-balance, requirement. */ 

V[TORQUE_N_WG] = v[TORQUE_K]; 

/* Calculate the cyclic torque. */ 
v[TORQUE_CYCLIC] = (cyclic_amplitude * 

Sin(((2.0 * v[POS_WGj) + cyclic_phase) * CONV)); 

/* Now apply the torque constraint imposed by the three-port gear-reduction 

* with position error included. Note that the torque on the flexspline 

* is defined to be positive in the negative rotation direction. */ 
v[TORQUE_N_FS] = -((-N / (1.0 - (N * v[DERF])) ) * v[TORQUE_N_WGj) ; 

V[TORQUE_N_CS] = ((N + 1.0) / (1.0 - (N * v[DERF]))) * v[TORQUE_N_WG]; 

/* Finally, calculate the torque on all of the harmonic-drive ports. */ 

V [ TORQ'JZ_WG ] = v [ TORQUE_K ] - v [ TORQUE_CYCLIC ] ; 
v[TORQUE_FS] = v[TORQUE_N_FS] - v[TORQUE_HD_FRICTION]; 
v[TORQUE_CS] = v[TORQUE_N_CS] - v[TORQUE_HD_FRICTION]; 

/* Finally, return the torques on the three hd ports. */ 

*torque_wg = v[T0RQUE_WG]; 

*torque_f s = v[TORQUE_FS] ; 

*torque_cs = v[TORQUE CS]; 


/* Function gear_tooth_hd_model() calculates the three-port harmonic-drive 

equations for the harmonic drive representation which uses inclined planes to 
* simulate the tooth-rubbing behavior inside the drive. */ 
void gear_tooth_hd_model(pos_wg, pos fs, pos cs, 

vel_wg, vel_fs, vel_cs, 
torque_wg, torque_fs, torque_cs) 
double pcs_wg, pos_fs, pos_cs; 
double vel_wg, vel_fs, vel_cs; 
double *torque_wg, *torque_fs, *torque_cs; 

{ 

—^_^^^£^_^J^_J^b,2^tb_ie_s__a_nd__po_si 1 tj.ons and the fs, cs, and wg ports 
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in the global data vector to begin calculations */ 
v[POS_WGJ = pos_wg; 
v[POS_FS] = pos_fs; 
v[POS_CS] = pos_cs; 
v[VEL_WG] = vel wg; 
v[VEL_FS] = vel_fs; 
v[VEL_CS] = vel_cs; 

/* F: '- rst ' calculate the continuity (or kinematic) constraints. 

* Find the veloc ity of the wave-generator with respect to the flexspline */ 
v[POS_WG_RELATIVE] = v[POS_WG] - v[POS FS]; 
v[VEL_WG_RELATIVE] = v[VEL_WG] - v[VEL _ FS]; 

/* Determine the kinematic error function and its derivative such that 
(d/dt)v[ERFN] = v[VEL_WG_RELATIVE] v[DERF] */ 
v(ERFN] = (N*tanl * ((error_amplitudeO * 

Sin(((1.0 * v[POS_WG_RELATIVE]) + error_phaseO) * CONV)) + 
(error_amplitudel * 

Sin (((2.0 * v[POS_WG_RELATIVE]) + error_phasel) * CONV)) + 

(error_amplitude 2 * 

Sin(((4.0 * v[POS_WG_RELATIVE]) + error_phase2) * CONV)))); 
v [DERF] - (N*tanl * ((1.0 * CONV * error_amplitudeO * 

Cos (((1.0 * v [ POS_WG_RELATIVE) ) + error__phaseO) * CONV)) + 
(2.0 * CONV * error_amplitudel * 

Cos(((2.0 * v[POS_WG_RELATIVE]) + error_phasel) * CONV)) + 
(4.0 * CONV * error_amplitude2 * 

Cos (((4.0 * v(POS_WG_RELATIVE)) + error_phase2) * CONV) ) )) ; 

/* Find the radial (vertical) movement of the tooth base. */ 
v[POS_TOOTH_BASE] = tanl * v[POS_WG_RELATIVE]; 
v [VEL_TOOTH_BASE] = tanl * v[VEL_WG_RELATIVE]; 

radial movement at the output of the position error element. */ 
v[POS_ERR] = v[POS_TOOTH_BASE] + V[ERFN J; 

v[VEL_ERR] = v[VEL_TOOTH_BASE] + (v[VEL_WG_RELATIVE] * v(DERF]); 

/* Find the radial movement of the tooth tip */ 
v[POS_TOOTH_TIP] = (v[POS_CS] - v[POS_FS]) /’tan2; 
v [VEL_TOOTH_TIP] = (v[VEL_CS] - v(VEL_FS]) / tan2; 

/ Find the velocity at the wave-generator surface. */ 
v[POS_WG_SURFACE] = ((cosl * v(POS WG]) + 

(sinl * v[POS_TOOTH_ BASE]) - 

(cosl * v[P0S_FS])); 

v[VEL_WG_SURFACE] = ((cosl * v(VEL WG]) + 

(sinl * v[VEL_TOOTH_BASE]) - 
(cosl * v[VEL_FS])); 

/ Find the velocity at the gear-tooth surface. */ 
v(POS_TOOTH_SURFACE] = ((cos2 * v(POS TOOTH TIP]) + 

(sin2 * v[POS_CS]) - 
(sin2 * v[POS_ FS])); 

v[VEL_TOOTH_SURFACE] = ((cos2 * v[VEL TOOTH TIP]) + 

(sin2 * v[VEL_CS]) - _ 

(sin2 * v[VEL FS])); 


/* Now ' calculate the Constitutive Relationships 

* Find the wave-generator friction. */ 
v[TORQUE_WG_FRICTION] = 

calculate_friction(v[VEL_WG_SURFACE], 

— b W C constant, b wql. b wg2. 
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0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 1000000 . 0 ); 

/* Find the friction at the tooth tip. */ 
v[TORQUE_TOOTH_TIP_FRICTION] = 

calculate_friction((v[VEL_TOOTH_BASE]- v[VEL_TOOTH_TIP]), 

b_tooth_tip_constant, b_tooth_tipl, b_tooth_tip2, 

0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 0 . 0 , 1000000 . 0 ); 

/* The total friction at the gear-tooth surface is due to the coulomb friction 

* {mu * N), the constant friction (b_tooth_surface_constant), the velocity- 

* dependent friction, and the cyclic friction. Calculate the component due 

* to everything except the coulomb friction. */ 

v[TORQUE_TOOTH_SURFACE_FRICTION] = 

calculate_friction(v[VEL_TOOTH_SURFACE], 

b_tooth_surface_constant, b_tooth_surfacel, b_tooth_surface2, 
stiction_torque_tooth, stiction_vel_tooth, 
v[POS_OUT], 

cyclic_friction_amp_tooth, cyclic_friction_phase_tooth, 

(max_friction_vel[model] [joint])) ; 


/* Find the harmonic drive stiffness coefficients: */ 
v[Kl] = (kl_constant + 

kl_amplitude * Sin(((2.0 * v (POS_WG__RELATIVE] ) + kl_phase) * CONV)); 

v[K2] = (k2_constant + 

k2_amplitude * Sin(((2.0 * v(POS_WG_RELATIVE]) + k2_phase) * CONV)); 
/* Calculate the torque on the spring. */ 

v [ TORQUE_K ] = ( < v [ K1 ] * (v[POS_ERR] - v [ POS__TOOTH_TI P ] ) ) + 

(v [K2 ] * Power ( (v [P0S__ERR] - v ( POS_TOOTH_TI P ) ) , 3))); 

/* Calculate the cyclic torque. */ 
v[TORQUE_CYCLIC] = (cyclic_amplitude * 

Sin (((2.0 * v[POS_WG_RELATIVE]) + cyclic_phase) * CONV) ) ; 

/* From the power conservation constraint on the position error element 

* calculate the torque seen on the input side. */ 

v[TORQUE_ERR_IN] = (1.0 + (v[DERFJ / tanl)) * v[TORQUE_K]; 

/* Set the sign of mu to control the direction of the coulomb friction force 

* along the gear-tooth surface. The sign of mu should be identical to the 

* sign of the surface rubbing-velocity. */ 
mu = copysign(mu_save, v[VEL_TOOTH_SURFACE]); 

/* If stiffness data is being collected, the friction acts in the opposite 

* direction as the joint is being loaded. */ 
if(stiffness_flag) 

mu = -mu_save; 

/* Now calculate the force/torque compatibility equations: 

* 

* Determine the output normal force. */ 

v[TORQUE_OUTPUT_NORMAL] = v[T0RQUE_K] + v[TORQUE_TOOTH TIP FRICTION]; 

/* For the sake of continuity in the simulation, ramp the constant friction 

* on the gear-tooth surface from zero to the desired value if at the 

* beginning of the step-response trial. This is a quick fix that should 

* be implemented in more general terms for all of the constant friction 

* components. */ 

if(!(stiffness_flag || stiffness_flag2)) 
if(current_time < 0.05) 

v[TORQUE_TOOTH_SURFACE_FRICTION] = ({((current_time - 0.05)/0.05) + 1.0) * 
_____ v[ TORQUE TOOTH SURFACE FRICTION]); 
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. . lnd the gear-tooth normal force. Note that this equation might need to 
be calculated iteratively if the velocity is zero since the coulomb • 
friction becomes a different function of the normal force. For now, I 
will start my simulations at a non-zero initial velocity so that the 
* velocity will never be zero. */ 

v[TORQUE_TOOTH_SURFACE_NORMAL] = {(v[TORQUE_OUTPUT_NORMAL] - 

(cos2 * v[TORQUE_TOOTH_SURFACE_FRICTION])) 
(sin2 + {mu * cos2))); 


/ 


± u- the gear-tooth surface-normal is negative then the mechanism through 
which the friction acts at the gear-tooth surface changes. In order to 
* describe this new mechanism, a different, but similar, model is needed 
Unfortunately, changing to the new model at this point in the simulation 
poses significant numerical problems. To avoid this mess, when the normal 
force becomes negative, I will set the coulomb friction to zero and 


* recalculate the tooth-surface normal force, 
if(v[TORQUE_TOOTH_SURFACE_NORMAL] < 0.0) 


'/ 


{ 


} 


mu = 0.0; 

v[TORQUE_TOOTH_SURFACE_NORMAL) = ((v[TORQUE_OUTPUT NORMAL) - 

(cos2 * v[TORQUE_TOOTH_SURFACE_FRICTION)) ) / 
(sin2 + (mu * cos2))); 


/* GlVe " th | s correct tooth-surface normal force, the coulomb friction and 
total friction at the tooth interface can be calculated. */ 
v[T0RQUE_T00TH_SURFACE_C0UL0M8] = mu * v(TORQUE_TOOTH SURFACE NORMAL); 
v[TORQUE_TOOTH_SURFACE_TOTAL_LOSS] = (v[TORQCJE_TOOTH_SURFACE_COULOMB] + 

-v[TORQUE_TOOTH_SURFACE_FRICTION]); 

/* Calculate the input normal force. */ 

v[TORQUE_INPUT_NORMAL] = v[TORQUE_ERR_IN] + v[TORQUE_TOOTH_TIP FRICTION); 

/* Now determine the normal force on the wave-generator. */ 
v[TORQUE_WG_NORMAL] = (1.0 / cosl) * (v[TORQUE_INPUT_NORMAL) + 

(v[TORQUE_WG_FRICTION] * sinl)); 

/ Find the normal force on the tooth tip. */ 

v [TORQUE_TOOTH__TIP_NORMAL ] = ((cos2 * v[TORQUE_TOOTH_SURFACE_NORMAL] ) - 

(sin2 * v(TORQUE_TOOTH_SURFACE_TOTAL LOSS])); 

/* Now the torques on the wave_generator, flexspline, and circular spline can 
* be determined. */ - 

V[TORQUE_WG] = ((v[TORQUE_CYCLICJ) + 

(v[TORQUE_WG_FRICTION] * cosl) + 

(v[TORQUE_WG_NORMAL] * sinl)); 
v [ T0RQUE_FS ] = ( (v[ TORQUE_TOOTH_TIP_NORMAL ]) - 
(v[TORQUE_WG_FRICTION] * cosl) - 
(v(TORQUE_WG_NORMAL) * sinl)); 
v[TORQUE_CS] = ((v[TORQUE_TOOTH_SURFACE_NORMAL] * cos2) - 

(v[TORQUE_TOOTH_SURFACE_TOTAL_LOSS] * sin2)); 

/ Finally, return the torques on the three hd ports. */ 

*torque_wg = v[TORQUE_WG]; 

*torque_f s = v[TORQUE_FS]; 

*torque_cs = v[TORQUE_CS]; 


'/***************-k****i,- k i'i' ici ' ici , i ' i ' itirici ' ic 

I/* 


'/* Functions which describe the energy of each h 
[/* 


*******************i r x** i ,i, i 'i, i , i ,* 1 , 1 , i '**i<**** 

armonic-drive model. 
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/*** ******************************************************************************** 

/ 

/* Function rotary_hd_energy() calculates the energy for the non-ideal, rotary 

* harmonic-drive model. It is called by the function odeint() after the completion 

* of each time step in order to integrate the individual energy components of the 

* harmonic drive model with tooth rubbing. This function takes the current time 

* step of the solver and integrates the energy over that time step using a simple 

* Euler integration scheme. As you might expect, these energy calculations can 

* sometimes have appreciable error, but they are useful, for the most part. */ 
void rotary_hd_energy(delta_t) 

double delta_t; 

{ 

double vel_wg_rad, vel_fs_rad, vel_cs_rad; 
double vel_hd_friction_rad; 
static double old_pos_k = 0.0; 

/* Be sure that all velocities and positions are in radians to ensure 

* proper unit-matching. */ 
vel_wg_rad = v[VEL_WG] * CONV; 
vel_fs_rad = v[VEL_FS] * CONV; 
vel_cs_rad = v[VEL_CS] * CONV; 

vel_hd_friction_rad = v[VEL_HD_FRICTION] * CONV; 

/* Calculate the potential energy stored in the spring by integrating 

* the spring force over the spring displacement. */ 

V[SPRING_ENERGY] = (v[SPRING_ENERGY] + 

((({v[POS_K] - v[POS_N_WG]) - old_pos_k) * CONV) * 
v[TORQUE_K]) ) ; - 

/* Calculate the cyclic torque energy by integrating the power over time. */ 
V[CYCLIC_ENERGY] = (v[CYCLIC_ENERGY] + 

(delta_t * vel_wg_rad * (-v [TORQUE_CYCLIC]))); 

/* Calculate the energy lost to friction by integrating the power dissipation 

* over time. */ 

V[HD_FRICTION_ENERGY] = (v[HD_FRICTION_ENERGY] + 

(delta_t * (v[TORQUE_HD_FRICTION] * vel_hd friction rad))); 

/* Calculate total potential energy in the system assuming the cyclic torque 

* can be treated like a spring. */ 

v[TOTAL_POTENTIAL_ENERGY] = v[SPRING_ENERGY] + v [CYCLIC_ENERGYj ; 

old_pos_k = (v[POS_K] - v[POS N WG]) ; 


/* Function gear_tooth_hd_energy() calculates the energy for the harmonic-drive 

* model which includes gear-tooth-rubbing effects. It is called by the function 
odeint() after the completion of each time step in order to integrate the 
individual energy components of the harmonic drive model with tooth rubbing. 
This function takes the current time step of the solver and integrates tne 
energy over that time step using a simple Euler integration scheme. As you 
might expect, these energy calculations can sometimes have appreciable error, 

* but they are useful, for the most part. */ 
void gear_tooth_hd_energy(delta_t) 

double delta_t; 

{ 

double vel_wg_rad, vel_fs_rad, vel cs rad; 
double vel_wg_surface_rad, vel_tooth_surface_rad; 
double vel_tooth_tip_rad, vel tooth base rad; 
static double old_pos_k = 0.0; 
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/* Be sure that all velocities and positions are in radians to ensure 

* proper unit-matching. */ 
vel_wg_rad = v(VEL_WG] * CONV; 
vel_fs_rad = v[VEL_FS] * CONV; 
vel_cs_rad = v[VEL_CS] * CONV; 
vel_wg_surface_rad = v[VEL_WG_SURFACE] * CONV; 
vel_tooth_surface_rad = v[VEL_TOOTH_SURFACE] * CONV; 
vel_tooth_tip_rad = v[VEL_TOOTH_TIP] * CONV; 
vel_tooth_base_rad = v(VEL_TOOTH_BASE] * CONV; 

/* Calculate the potential energy stored in the spring by integrating 
the spring force over the spring displacement. */ 
v[SPRING_ENERGY] = (v[SPRING_ENERGY] + 

( ( ( (v [POS_ERR] - v[ POS_TOOTH_TIP ]) - old_pos_k) * CONV) * 
v[TORQUE_K])) ; 

/* Calculate the cyclic torque energy by integrating the power over time. */ 
v[CYCLIC_ENERGY] = (v[CYCLIC_ENERGY] + 

(delta_t * vel_wg_rad * (-v[TORQUE_CYCLICj))); 

/* Calculate the energy lost to friction by integrating the power dissipation 

* over time. */ 

v[WG_FRICTION_ENERGY] = (v[WG_FRICTION_ENERGY] + 

(delta_t * (v[TORQUE_WG_FRICTION] * vel wg surface rad))) 
v[TOOTH_TIP_FRICTION_ENERGY] = (v(TOOTH_TIP_FRICTION_ENERGY] +~ 

(delta_t * 

(v[TORQUE_TOOTH_TIP_FRICTION] * 

(vel_tooth_base_rad - vel tooth tip rad)))); 
v[TOOTH_SURFACE_FRICTION_ENERGY] = (v[TOOTH_SURFACE_FRICTION_ENERGY] +■ 

<delta_t * 

(v[TORQUE_TOOTH_SURFACE_FRICTION] * 

vel_tooth_surface_rad) ) ) ; 

v[TOOTH_SURFACE_COULOMB_ENERGY] = (v(TOOTH_SURFACE_COULOMB ENERGY] + 

(delta_t * 

(V[TORQUE_TOOTH_SURFACE_COULOMB] * 
vel_tooth_surface_rad))); 

v[TOOTH_SURFACE_TOTAL_LOSS_ENERGY] = (V[TOOTH_SURFACE_COULOMB_ENERGY) + 

v[TOOTH_SURFACE_FRICTION ENERGY]); 

/* Calculate total potential energy in the system assuming the cyclic torque 
* can be treated like a spring. */ 

v[TOTAL_POTENTIAL_ENERGY] = V[SPRING_EN£RGY) + v(CYCLIC ENERGY); 

°id_pos_k = (v[POS_ERR) - v[POS TOOTH TIP]); 


/*■************* 

/* 

/* Functions 
/* 


****■************************************■*■***********. 




which calculate model parameters based on input values. 


****************** 






/* Function calculate_rotary_params() is used to take the 
values and determine relevant parameters for the simple 
* harmonic-drive model. */ 
void calculate_rotary_params() 


input parameter 
non-ideal rotary 


double reduction_factor; 
double frequency; 
double tempi, temp2; 
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/* Reflect the non-linear stiffness coefficients as measured from the output 

* side of the transmission to the input side. */ 
kl_constant = (kl_output_constant / Power(N, 2)); 
kl_amplitude = (kl_output_amplitude / Power(N, 2)); 
kl_phase = kl_output_phase; 

k2_constant = (k2_output_constant / Power(N, 4)); 
k2_amplitude = (k2_output_amplitude / Power (N, 4)); 
k2_phase = k2_output_phase; 

/* Calculate the coulomb friction in the harmonic drive based on the width 
of the hysteresis loop in the harmonic-drive stiffness profile. */ 
b_hd_constant = (hysteresis_width / 2.0); 

/* Calculate the reduction ratio between the input rotation and the velocity 
seen by the friction element mounted between the flexspline and circular 

* spline ports. */ 
if(joint == WRIST) 

reduction_factor = N; 
else 

reduction_factor = (N + 1.0); 

/* Calculate the dynamic friction coefficients by subtracting the input 
and output damping components from the total dynamic friction using 
conservation of power. These coefficients should be positive. */ 
b_hdl = ((Power(reduction_factor, 2)) * (b_totall - b_in)) - b_out; 
b_hd2 = ((Power(reduction_factor, 4))) * b_total2; 

/ Convert the stiction and cyclic friction from input rotation units to 
output rotation. Don't upset the signs of these values. */ 
stiction_torque_hd = reduction_factor * stiction torque; 
stiction_vel_hd = (stiction_vel / reduction factor); 

cyclic_friction_amp_hd = reduction_factor * cyclic friction amplitude; 
cyclic_friction_phase_hd = cyclic_friction phase; 

/* Calculate the natural frequency of the non-rigid vibrational mode of the 

* two-mass system assuming linear stiffness. */ 
tempi = ((kl_constant / CONV) * 

((1.0 /inertia_in) + 

(1.0 / (inertia_out / (Power(reduction_factor, 2.0)))))); 
temp2 = (1.0 / (Power(2.0, 1.5) * PI)); 
frequency = (temp2 * Power((2.0 * tempi), 0.5)); 

printf("\n\n The natural frequency of the non-rigid vibrational mode\n"); 
printf ( for a two—mass linear system is; %f Hz,\n", frequency); 
printf( Therefore, resonance vibration should appear at input rotational\n"); 
printf( velocities of: %f input_deg/sec.\n", ((frequency * 360.0) / 4.0)); 

printf(" %f input_deg/sec.\n", ((frequency * 360.0) / 2.0)); 

printf(" %f input_deg/sec.\n", ((frequency * 360.0) / 1.0)); 


/ Print out these new parameters if in debugging mode. */ 
if(debug) 


printf("\n b_hd_constant : 
printf("\n b_hdl: 
printf("\n b_hd2: 
printf("\n kl_constant: 
printf("\n k2_constant: 
printf("\n") ; 


%4.201f", 
%4.20 If ", 
%4 . 201f ", 
%4.201f", 


b_hd_constant) ; 
b_hdl); 
b_hd2); 
k1 constant); 


%4.201f\n", k2 constant); 


/ Function calculate_gear_tooth_params() is used to take the input parameter 
values and determine: (1) trigonometric and geometric relationships. (2) the 
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* tooth surface friction coefficients, and (3) the ideal spring stiffness 

* coefficients. */ 

void calculate_gear_tooth_params(’ 

{ 

double stiffness_factor_A; 
double stiffness_factor B; 
double wg_to_tooth_factor; 
double damping_factor A; 
double damping_factor B; 
double damping_factor C, 
double damping_factor D; 
double inertia_in_eff, inertia_out_eff; 
double tempi, temp2, frequency; 

/* Initialize the coulomb friction coefficient. */ 
mu = mu_save; 

/* Calculate the gear-tooth trig values. */ 
tan2 = tan(CONV * tooth angle); 
sin2 = sin(CONV * tooth angle); 
cos2 = cos(CONV * tooth angle); 

/* The wg_angle can be found from the gear ratio and the tooth angle 
* using the relationship: (1/(N+l)) = tanl tan2. */ 
wg_angle = atan(1.0 / (tan2 * (N + 1.0))); 
tanl = tan(wg_angle); 
sinl = sin(wg_angle); 
co s1 = cos(wg_angle); 

/* ”° W Calculate the stiffness of the spring that deforms vertically 
from the rotational stiffness provided. */ 
stiffness_factor_A = (cos2 + (mu * sin2)) / (sin2 - (mu * cos2)); 
stiffness_factor_B = -(1.0 / (tanl - (1.0 / tan2))); 
kl_constant = (kl_output_constant * 

(stiffness_factor_B / (stiffness_factor A - tanl) ) ) ; 
kl_amplitude = (kl_output_amplitude * 

(stiffness_factor_B / (stiffness_factor A - tanl))); 

^l_phase = kl_output_phase; 
k2_constant = (k2_output_constant * 

(Power(stiffness_factor_B, 3) / (stiffness factor A - tanl))); 
k2_amplitude = (k2_output_amplitude * 

(Power(stiffness_factor_B, 3) / (stiffness factor A - tanl))); 
k2_phase — k2_output phase; 

/* Fr ° m thS h V steres i s width on the stiffness curve and the values of the 
other constant friction components in the harmonic drive, determine 

e appropriate constant friction component on the gear tooth surface */ 
b_tooth_surface_constant = ((1.0 / ((stiffness_factorA * cos2) + sin2)) * 

((hysteresis_width / 2 . 0 ) - 

((stiffness_factor_A - tanl) * b_tooth_tip constant) 
(((sinl * tanl) + cosl) * b_wg constant)) ) ; 

/* Scale the stiction torque and cyclic friction torque from the input 
side (wg) to the gear tooth surface. */ 

W9 T t0 _t° 0th - factor = (1.0 / ((tanl * cos2) + (tanl * sin2 * tan2))); 
stiction_torque_tooth = wg_to_tooth_factor * stiction torque; 
stiction_vel_tooth = (stiction_vel / wg_to_tooth_factor); 

^t 1C -f^r i0n - a : P - t0 ° th = w< 3_ to _tooth_factor * cyclic_friction_ampiitude; 
cyclic_fnction_phase_tooth = cyclic_f riction_phase; 

/* By applying conservation of power to the damping components in the 
model, the linear and cubic coefficients for the gear-tooth-surface 

* ? 3n bS f ° Und by subtrac ting the known damping components from 
-t he total joint dampingjnd matching coefficients. This has to be 
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* done differently for configuration 1 and configuration 2 since the 

* dampers outside the harmonic drive depend on different velocities. */ 
damping_factor_A = ( (cos2 / tan2) + sin2)/ 

damping_factor_B = (cosl + (sinl * tanl)); 
damping_factor_C = (1.0 / (N + 1.0)); 
damping_factor_D = (1.0 + (1.0 / N)); 
if(joint == WRIST) 


b_tooth_surfacel = ((b_totall 

b in - 


b tooth surface2 


(b_wgl * (Power((damping_factor_B * damping_factor_D), 2))) 
(b_out * (Power((-1.0 / N), 2)))) / 

(Power((damping_factor_A * (1.0 / N) ) , 2))); 

= ( (b_total2 - 

(b_wg2 * (Power((damping_factor_B * damping_factor D), 4))) 
(Power((damping_factor_A * (1.0 / N) ) , 4))); 


else 

{ 


b_tooth_surfacel = ((b_totall - 

b in - 


(b_wgl * (Power((damping_factor_B) 
(b_out * (Power((damping_factor C) 


2 ) ) ) - 
2 ) ) ) ) / 


(Power((damping_factor_A * damping_factor_C), 2))); 
b_tooth_surface2 = ((b_total2 - 

(b_wg2 * (Power((damping_factor_B) , 4)))) / 

(Power((damping_factor_A * damping_factor_C) , 


4 ) ) ) ; 


/* Calculate the natural frequency of the non-rigid vibrational mode 

* two-mass system assuming linear stiffness. First calculate the ef 
input and output inertias seen by the internal spring, then use th 

* ideal linear stiffness value, with friction removed, to determine 

* resulting natural frequency of the linear system. */ 
inertia_in_eff = (inertia_in / (tanl * tanl)); 
inertia_out_eff = (inertia_out * tan2 * tan2); 

tempi = ((kl_constant / CONV) * 

((1.0 /inertia_in_eff) + (1.0 / inertia_out_eff))); 

temp2 = (1.0 / (Power(2.0, 1.5) * PI)); 
frequency = (temp2 * Power((2.0 * tempi), 0.5)); 

printf("\n\n The natural frequency of the non-rigid vibrational mode\ 
print f ( for a two-mas s linear system is: %f Hz.\n", frequency); 
printf ( Therefore, resonance vibration should appear at input rotati 
printf (" velocities of: %f input_deg/sec.\n", ((frequency * 360.0) / 

print f ( " %f input_deg/sec.\n", ((frequency * 360.0) / 

print f ( " %f input_deg/sec.\n", ((frequency * 360.0) / 

/* Print out these new parameters if in debugging mode. */ 
if(debug) 


of the 
fective 


ona 1\n " ) ; 
4.0)); 
2 . 0 )); 
1 . 0 )); 


printf("\n wg_angle: %4.201f deg", (wg_angle / CONV)); 

printf( \n b_tooth_surface_constant: %4.201f", b_tooth_surface constant); 
printf( \n b_tooth_surfacel: %4.201f", b tooth surfacel); 

printf("\n b_tooth_surface2: %4.201f", b2tooth~surface2); 

printf("\n kl_constant: %4.201f", kl_constant); 

printf( \n k2_constant: %4.201f\n", k2 constant); 

printf("\n"); 

printf("\n tooth surface to input scaling factor: %4.201f\n", 

(1.0 / (damping_factor_A * damping_factor_C))); 
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z^****************************************************.****^*^ 

7 * 

/* Functions which model the friction and amplifiers and collect stiff 
/ * 

/*********************-k*1,1c**-kit-k1,**ir1c-k***-l,-iri,i,^i,i r i r - k i,i ri(i(iritiri , 


**-k'kic*-ki(icic-Xic*ic*-k-x*ic-k* 


w / 

ness data. */ 


'/ 


*★★*★★■****■******.* 


*********: 


/* Function irequlateO takes as arguments the present amp current, the requested 

amp current, and the current motor velocity and returns the actual current that 
* the amps can produce. */ 

double iregulate(ipresent, irequested, omega) 
double ipresent; 
double irequested; 
double omega; 

{ 

double voltage_upper_limit; 
double voltage_lower_limit; 
double current_upper_limit; 
double current_lower_limit; 
double currentl, current2; 
double ireturn; 

/* Make sure that the current requested is within the amp's operating range */ 
if (irequested > IMAX) 

( 

fprintf(stderr, "\n WARNING: requested current is above allowed amp range'\n")■ 
printf("\n Current Time: "); 
irequested = IMAX; 

} 

else if (irequested < -IMAX) 

( 

fprintf(stderr, "\n WARNING: requested current is below allowed amp ranqe'Xn")- 
printf("\n Current Time: "); 
irequested = -IMAX; 

) 

/* At the present operating current, ipresent, calculate the maximum 
allowable voltage that the amps can provide as calculated from 
the voltage-current curves for the amps in the Aerotech catalog. 

Assume that, if the current is negative, the maximum voltage that me 
amp can reach is the saturation voltage, VMAX. Also assume that the 
minimum voltage for negative current is described by a curve similar 
to the voltage-current for positive current and voltage. Lastly, tne 
minimum voltage for positive current is equal to the negative voltage 
* threshold of the amp, -VMAX. */ 
if(ipresent >= 0) 


{ 


) 

else 

( 


if (ipresent <= IMAX) 

voltage_upper_limit = VMAX + VI_SLOPE * ipresent; 
else 

voltage_upper_limit = VMAX + VI_SLOPE * IMAX; 
voltage_lower_limit = -VMAX; 


if (ipresent >= -IMAX) 

voltage_lower_limit = -VMAX + VI_SLOPE * ipresent; 
else 

voltage_lower_limit = -VMAX + VI_SLOPE * (-IMAX); 
voltage_upper_limit = VMAX; 


—. Now from this maximum or m inimum allowed voltage at present current, caiculate 
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* the allowable current range that the requested current can achieve by assuming 

* that the amp voltage is equal to the voltage across the motor (= Kb * velocity 

* or back EMF) plus the voltage drop due to current flowing through the' amp 

* resistance. */ 

current 1 = (voltage_upper_limit - (motor_kb * omega))/amp resistance; 
current2 = (voltage_lower_limit - (motor_kb * omega))/amp_resistance; 
if(currentl < current2) 

{ 

current_upper_limit = current2; 
current_10wer_limit = currentl; 


else 


current_upper_limit = currentl; 
current_lower_limit = current 2 ; 


/* Now, if the requested current, irequested, is between the upper and lower 

* current bounds, return the value, otherwise return the upper or lower 

* limit. */ 

if((irequested <= current_upper_limit) && (irequested >= current lower limit)) 
ireturn = irequested; 

else 

{ 

if((fabs(current_lower_limit - irequested)) < 

(fabs(current_upper_limit - irequested))) 
ireturn = current_lower_limit; 
else 

ireturn = current_upper_limit; 

} 

/* Now add a relative damping factor to the current. Make it behave so that 
the current can't change from the old to the new value instantaneously, bur, 
instead, a part of the new value, ireturn, is averaged with part of the old 
value, ipresent, to create the actual new current somewhere between the two 
values. This step is necessary to keep the amps current from fluctuating 
rapidly between each subsequent step and possibly going unstable. The AMP 
DAMPING factor is a number between 0 and 1.0 that represents the amount 

* of damping in the amps current response. Note that since step size can vary 
significantly, and that this damping factor is calculated on every step, the 
amount of damping will vary with the changing step size. */ 

ireturn = (AMP_DAMPING * ipresent) + ((1.0 - AMP_DAMPING) * ireturn); 

/* Give a run-time warning if the current exceeds the amp's allowed range. */ 

if ((ireturn > IMAX) || (ireturn < -IMAX)) 

{ 

fprintf(stderr, "\n WARNING: allowable amp current range exceeded!\n"); 
print f ( "\n Current Time: "); 


return(ireturn) ; 


Function calculate_friction() takes 10 arguments: (1) the velocity at the 
friction interface, (2,3, and 4) the constant, linear, and cubic veiocity- 
dependent friction coefficients, (5) the stiction torque, (6) the velocity 
where the stiction is deactivated, (7) the reference position for the cyclic 
friction, (8 and 9) the amplitude and phase of a sinusoidal friction function, 
and (10) the velocity for which the damping curve begins to decrease. 

This function returns the resulting value of the friction given these 
parameters. Note that the friction amplitude sh ould not be larger than 
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b_constant or else the resulting friction can become negative. Additionally, 

* if the velocity is very small, this friction model is undefined, and a warning 
is given. Finally, if the velocity is too large, the cubic term of the 

velocity-dependent damping may decrease and become negative. To prevent this 

* the velocity dependent damping is calculated for omega_max at all velocities 

* greater than omega max. */ 

|double calculate_friction(omega, 

b_constant, bl, b2, 
stiction_torque, stiction_vel, 
theta, 

friction_amplitude, f riction_phase, 
omega_max) 

double omega; 

double b_constant, bl, b2; 
double stiction_torque, stiction_vel; 
double theta; ~ 

double friction_amplitude, friction phase; 
double omega_max; 

double friction, damping; 


/ Make sure that the velocity is not near zero. */ 
if (((fabs(omega)) < TINY_VEL) && 

( (stiction_torque + friction_amplitude + b_constant) 
(!stiffness_flag) && (!stiffness_flag2)) 


= 0.0) i& 


{ 

fprintf(stderr, 
fprintf(stderr, 
|undefined!\n") ; 

return(0.0); 


"\n ERROR: Function calculate_friction() .") ; 

"\ n Velocity is near zero and the model is 


/* If the stiffness data is being collected and the joint is being loaded, 

* K ^ etUrn the ve l° c ity - independent friction as if the velocity were positive. 
Note: do not use the stiction torque. */ 
if (stiffness_flag) 


{ 


friction = (-b_constant + 

(-friction_amplitude * 

Sin ( (theta + friction_phase) * CONV))); 
return(friction); 


/* If stiffness c * ata is being collected and the joint is being unloaded, 

return the velocity-independent friction as if the velocity were positive */ 
if (stiffness_flag 2 ) 


{ 


friction = (b_constant + 

(friction_amplitude * 

Sin ((theta + friction_phase) * CONV))); 
return(friction); 


} 


/* Make sure that the cylic friction amplitude is less than the constant friction 
* coefficient. */ 

(friction_amplitude > b constant) 

{ 

fprintf(stderr, "\n ERROR: Function calculate_friction().") ; 

} f P rintf (stderr, "\n Cyclic friction amplitude too large!\n") ; 


/* Calculate the static and constant friction as follows: 

* lf t [j e y elocit y is below the stiction_vel, calculate the friction oaseo 
- the llnear flt between max imum stiction at zero velocity and b constant 


218 



Appendix B: S imulalion Code 


* at stiction_vel. Otherwise the friction is equal to b_constant in the 

* direction opposing the velocity. Add the cylic friction accordingly. */ 
if (omega >= stiction_vel) 

friction = (b_constant + 

(friction_amplitude * 

Sin((theta + friction_phase) * CONV))); 
else if ((omega > 0.0) && (omega < stiction vel)) 
friction = (stiction_torque + 

(( stiction_torque - b_constant)/(-stiction_vel) ) * omega); 
else if (omega <= (-stiction_vel)) 
friction = ((-b_constant) + 

(-friction_amplitude * 

Sin((theta + friction_phase) * CONV))); 

else 

friction = ((-stiction_torque) + 

(((~b_constant) - (-stiction_torque))/(-stiction_vel)) * omega); 

/* Now calculate the velocity-dependent friction. */ 
if(fabs(omega) > omega max) 

damping = ((bl * copysign(omega_max, omega)) + 

(b2 * Power((copysign(omega_max, omega)), 3))); 

else 

damping = (bl * omega) + (b2 * Power(omega, 3)); 

/* Return the aggregate friction value. */ 
return (friction + damping); 


/* Function collect_stiffness_data() can be used to collect a set of output 

* position versus output torque data. This function uses the specified 
harmonic drive function to calculate the resulting torques when a displacement 

* is applied to the joint output. The joint is configured with the wave 
generator locked to the circular spline. While the joint is being loaded, 

* data is collected in position increments of pos_increment until the maximum 
torque is reached. At that point, the harmonic drive model is switched since 

* the friction changes direction and the remaining data is collected. A positive 
stiffness curve is returned for positive displacement. Assuming directional 

* symmetry in the drive, this curve can be flipped about the x- and y- axes to 
determine the shape of curve for negative displacements. Be sure to load the 
drive in the direction that ensures that the gear-tooth-surface normal force 

* remains positive. */ 
void collect_stiffness_data() 

{ 

int i; 

double torque[STIFFNESS_POINTS]; 

double displacement(STIFFNESS_POINTS]; 

double pos_increment; 

double actual_disp, torque_wg, torque_fs, torque cs; 

char unix_command[200); 

int points = STIFFNESS_POINTS; 

static double max_disp[NUM_OF_AXIS) = (0.1, 0.4, 0.35); 

static char *stiffness_filename[NUM_OF_AXIS] = {"shld_stiffness data.dat", 

"elb_stiffness_data.dat", 
"wrist_stiffness_aata.dat"}; 

FILE *outfile, *gnufile; 


/* Calculate the position increment. */ 

P os _i cerement = max_disp[joint] / ((double) (points/2.0)); 

/* Begin loading the output link of the joint in compression. */ 
for(i = 1; j <= (points/2); i++) 
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actual_disp — ((double) <i— 1 )) * pos increment; 

if ( joint == WRIST) 

{ 

(*(hd_model[model])) ( 0 . 0 , actual_disp, 0 . 0 , 0 . 0 , 0.0, 0 . 0 , 

&torque_wg, &torque_fs, Storque_cs); 
torque[i] = torque fs; 

} 

else 

{ 

actual_disp = -actual_disp; 

(*(hd_model[model]))(actual_disp, 0.0, actual_disp, 0.0, 0.0, 0.0, 

&torque_wg, &torque_fs, Storque cs); 
torque[i] = torque_cs; 

} 

displacement[i] = fabs(actual_disp) ; 
if(debug) 

print_variables() ; 


/* Adjust the flags to allow proper friction values to be computed. */ 
stiffness_flag = FALSE; 
stiffness_flag 2 = TRUE; 

/* Collect data while unloading. Note the appropriate harmonic drive 
* model for the reverse-friction case is selected. */ 
for(i = ((points/ 2 ) + 1 ); i <= points; i++) 

{ 

actual_disp = ((double) (points - i)) * postincrement; 

if ( joint == WRIST) 

{ 

(*(hd_model[model])) ( 0 .0, actual_disp, 0.0, 0.0, 0 . 0 , 0.0, 

&torque_wg, &torque__fs, &torque_cs) ; 
torque[i] = torque_fs; 

} 

else 

{ 

actual_disp = -actual^disp; 

(*(hd_model[model]))(actual_disp, 0.0, actual_disp, 0.0, 0.0, 0.0, 

&torque_wg, 4 torque_fs, &torque_cs); 
torque[i] = torque_cs; 

) 

displacement[i] = fabs(actual_disp); 
if(debug) 

print_variables(); 

} 


/* Now store and print the data. */ 
create_file(&outfile, outdatapath, stiffness 
export__vectors (outfile, torque, displacement, 
close_file(out file); 


filename[joint]) 
points, stiffne 


ss 


filename[joint 


if(fast_plot) 

{ 

create_file(&gnufile, outdatapath, gnuplot_filename) ; 
fprintf(gnufile, "cd '%s'\n", outdatapath); 
fprintf(gnufile, "set xlabel 'Torque (N*m)'\n"); 
fprintf(gnufile, "set ylabel 'Displacement (deg)'\n"); 
fprintf(gnufile, "set title 'Stiffness Curve'\n"); 
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fprintf(gnufile, "plot , %s , \n", stiffness_filename[joint]); 
fprintf(gnufile, "pause -l\n"); 
close_file(gnufile); 

sprintf(unix_command, "gnuplot %s%s", outdatapath, gnuplot_filename); 
printf("\n Executing the unix command: %s\n", unix command); 
system(unix command); 


else 

{ 


sprintf(unix_command, "xgraph %s%s -t Stiffness_Curve\ 

-y Displacement_deg -x Torque_N*m i", 
outdatapath, stiffness_filename[joint]); 
printf("\n Executing the unix command: %s\n", unix command); 
system(unix command); 


exit ( 0 ); 


/******************************* *********************■»****■,*»**,**»,,,******, 

/* Functions which input and output data. «/ 

/ * 

' / 

/*****************************************.***********, ************,***,***,*.**,,, 

/* Function get_arguments() parses the command-line arguments. */ 
void get_arguments(int argnum, char **args) 

{ 

int arg_index; 

for (arg_index = 1 ; arg_index < argnum; arg_index++) 

{ 

if (!strcmp("-debug", args[arg^index])) 
debug = TRUE; 

else if (!strcmp("-energy", args[argindex])) 

{ 

print_energy = TRUE; 
calc_energy = TRUE; 

) 

else if (!strcmp("-stiffness", args[arg_index)) ) 
stiffness_flag = TRUE; 

else if (!strcmp("-ideal", args(arg_index] ) ) 
ideal_flag = TRUE; 

else if (!strcmp("-plot", args[arg_index])) 
plot = TRUE; 

else if (!strcmp("-h", args(arg_index)) || 

!strcmp("-?", args[arg_index]) || 

!strcmp("-help", args(arg_index))) 

{ 

fprintf(stderr, "\n Program %s: ",args[ 0 ]); 
usage() ; 
exit( 0 ); 

) 

else 

{ 


args[ 0 ]) ; 


fprintf(stderr,"\n ERROR while executing %s: arguments are wrong!\n", 


usage (); 
exit ( 1 ) ; 
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/* Function usage() prints out the proper usage of the command-line options. */ 
void usage() 


fprintf(stderr, 
fprintf(stderr, 
fprintf(stderr, 
fprintf(stderr, 
fprintf(stderr, 
fprintf(stderr, 
currents\n"); 
fprintf(stderr. 


\n Valid command-line options:\n"); 


-debug : 
-energy: 
-stiffness 
-ideal: 
-plot: 


activates debugging mode\n"); 
prints energy during runtime\n"); 
collects and plots solely stiffness data\n 
imposes an ideal haromonic-drive model\n") 
save plotting data for many input 


-help : prints out help messageXn"); 


/* Function read_input_data() parses all of the values in the input data file. */ 
void read_input_data(file_pointer) 

FILE *file pointer; 


char temp_name[50); 
double temp_number; 

while(fscanf(file_pointer, "%s %lg %*s\n", temp_name, &temp_number) != EOF) 

if (0 == strcmp("inertia_in", temp_name)) 
inertia_in = temp_number; 
else if (0 == strcmp("inertia_out", tempjame) ) 
inertia_out = temp_number; 

else if (0 == strcmp("kl_output_constant", temp_name)) 

^l_output_constant = temp number; 
else if (0 == strcmp("kl_output_amplitude", temp_name)) 
kl__output_amplitude = temp number; 
else if (0 == strcmp("kl_output_phase", temp_name)) 
kl_output_phase = temp_number; 
else if (0 == strcmp("k2_output_constant", temp name)) 
k2_output_constant = temp number; 
else if (0 == strcmp("k2_output_amplitude", temp name)) 
k2_output_amplitude = temp number; 
else if (0 == strcmp("k2_output_phase", terap_name)) 
k2_output_phase = temp number; 
else if (0 == strcmp("hysteresis_width", temp_name)) 
hysteresis_width = temp_number; 
else if (0 == strcmp("b_totall", temp_name)) 
b_totall = temp_number; 
else if (0 == strcmp("b_total2", temp_name)) 
b_total2 = temp_number; 
else if (0 == strcmp("b_in", temp_name)) 
b_in = temp__n umber; 

else if (0 == strcmp("b_wg_constant", temp_name)) 
b_wg__constant = temp number; 
else if (0 == strcmp("b_wg1", temp name)) 
b_wgl = temp_number; 

else if (0 == strcmp("b_wg2", temp name)) 
b_wg2 = temp_number; 

else if (0 == strcmp("b_tooth_tip_constant", temp name)) 
b_tooth_tip_constant = temp number; 
else if (0 == strcmp("b_tooth_tipl", temp_name)) 
b_tooth_tipl = temp number; 
else if (0 == strcmp("b_tooth_tip2", temp_name)) 

_b__t_ooth tip2 = temp number; 
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else if (0 == strcmp("cyclic_friction_amplitude", temp_name)) 
cyclic_friction_amplitude = temp_number; 
else if (0 == strcmp("cyclic_friction_phase", temp_name)) 
cyclic_friction_phase = temp_number; 
else if (0 == strcmp("stiction_torque“, temp_name)) 
stiction_torque = temp_number; 
else if (0 == strcmp("stiction_vel", temp_name)) 
stiction_vel = temp_number; 
else if (0 == strcmp("b_out", temp_name)) 
b_out = temp_number; 
else if (0 == strcmp("mu", temp_name)) 
mu_save = temp_number; 

else if (0 == strcmp("tooth_angle" / temp_name)) 
tooth_angle = temp_number; 
else if (0 == strcmpC'N", temp_name) ) 

N = temp_number ; 

else if (0 == strcmp("cyclic_amplitude", temp_name)) 
cyclic_amplitude = temp_number; 
else if (0 == strcmp{"cyclic_phase" / temp_name)) 
cyclic_phase = temp_number; 
else if (0 == strcmp("error_amplicudeO", temp_name)) 
error_amplitudeO = temp number; 
else if (0 == strcmp("error_phaseO", temp_name)) 
error_phaseO = temp_number; 
else if (0 == strcmp("error_amplitudel", temp_name)) 
error_amplitudel = temp__number ; 
else if (0 == strcmp("error_phasel", temp_name)) 
error_phasel = temp_number; 
else if (0 == strcmp{"error_amplitude2", temp^narae)) 
error_amplitude2 = temp_number; 
else if (0 == strcmp{"error_phase2", temp_name)) 
error_phase2 = temp_number; 
else if (0 == strcmp("motor_kt", temp_name)) 
motor_kt = temp_number; 
else if (0 == strcmp("motor_kb", temp_name)) 
motor_kb = temp_number; 

else if (0 == strcmp("amp_resistance", temp_name)) 
amp_resistance = temp_number; 
else if (0 == strcmp("irequested", temp_name)) 
irequested = temp_number; 

/* Initial conditions, */ 

else if (0 == strcmp("initial_pos_in", temp_name)) 
xstart_sav[POS_IN_INDEX] = temp_number; 
else if (0 == strcmp("initial_pos_out", temp_name)) 
xstart_sav[POS_OUT_INDEX] = temp number; 
else if (0 == strcmp("initial_vel_in", temp_name)) 
xstart_sav(VEL_IN_INDEX] = temp number; 
else if (0 == strcmp("initial_vel_out", temp_name)) 
xstart_sav[VEL_OUT_INDEX] = temp number; 

/* Time conditions */ 

else if (0 == strcmp("initial_time“, temp_name)) 
initial_time = temp_number; 
else if (0 == strcmp("final_time", temp_name)) 
final_time = temp_number; 
else if (0 == strcmp("step", temp_name)) 
step = temp_number; 

/* Fast plotting flag. */ 

else if (0 == strcmp("fast_plot", temp_name)) 
fast_plot = ((int) temp number); 


/* Desi red harmonic-drive model. */ 
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else if (0 == strcmp("model", temp_name)) 
model = ((int) temp number); 


else 

{ 


fprintf(stderr,"\n Couldn't parse parameter %s\n", temp name); 


/* Function read_input_filenames() reads the data from the second input file 

* as First it checks to see if the first value on each line is 1. 

If the va l ue is 1/ then it stores the line number in the vector index) j 

* so that it can be used later to remember which output variables to store. 
Additionally, this function also stores the other informaiton from each line 
that has a 1 value to be used later for graph titles and output filenames. 

As this function steps through each line of the input file, it counts the 
number of output variables selected and stores the result in the variable 

* num_selected_output_var. */ 

void read_input_filenames(file pointer) 

FILE *file_pointer; 

{ 

int line_number = 0; 
int temp_flag; 
char temp_title [50); 
char temp_units(50]; 

while(fscanf(file_pointer, "%d %s %s\n", &temp_flag, temp_title, temp^units) != EOF) 

if (temp_flag) 

{ 

index[num_selected_output_var] = line_number; 
strcpy(graph_title[num_selected_output_var] , temp title); 
strcpy(graph_units[num_selected_output_var], temp_units); 
sprintf(output_datafile[num_selected_output var], "%s hd %s.dat", 
joint_name[joint], temp_title); 
num_selected_output_var++; 


line__n umber + +; 


/* Function print_input_parameters() prints out the names and values of 
* the input parameters. */ 
void print_input_parameters() 


print f ("\n"); 

printfC Model Parameters :\n") ; 
printf(" Input Inertia 
printfC Output Inertia 
printfC kl constant 
printfC kl amplitude 
printfC kl phase 
printfC k2 constant 
printfC k2 amplitude 
printfC k2 phase 
printfC Hysteresis Width 
printfC b_totall 
printf(" b total2 


%4 . 101f\n", 
%4 . lOlf\n", 
%4 . 101f\n", 
%4 . lOlf\n", 
%4 .101f\n", 
%4 . 101f\n", 
%4 .101f\n", 
%4.101f\n", 
%4.101f\n", 
%4 . 101f\n", 
%4 .101f\n" 


inertia_in); 
inertia_out); 
kl_output_constant); 
kl_output_amplitude); 
kl_output_phase); 
k2_output_constant); 
k2_output_amplitude); 
k2_output_phase); 
hysteresis_width); 
b_totall) ; 
b total2) 
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printf (" b_in 
printf (" b_wg_constant 
printf(" b_wgl 
printf (" b_wg2 

printf(" b_tooth_tip_constant 
printf(" b_tooth_tipl 
printf(" b_tooth_tip2 
printf{" Cyclic Friction Amp 
printf(" Cyclic Friction Phase = 
printf(" Stiction Torque 
printf(" Stiction Vel 
printf(" b_out 
printf(" mu 
printf(" tooth_angle 
printf{" Gear Ratio, N 
printf(" Cyclic Amplitdue 
printf(" Cyclic Phase 
printf(" Error Amplitude 0 
printf(" Error Amplitude 1 
printf (•' Error Amplitude 2 
printf(" Error Phase ShiftO 
printf(" Error Phase Shiftl 
printf(" Error Phase Shift2 
printf(" Motor Torque Kt 
printf(" Motor Back EMF Kb 
printf(" Motor Resistance 
printf(" Motor Current 
printf ("\n") ; 

printf("Initial Conditions:\n"); 
printf(" Input Position 
printf(" Output Position 
printf(" Input Velocity 
printf(" Output Velocity 
printf("\n"); 
printf(" Initial Time 
printf(" Final Time 
printf(" Step Size 


%4 .lOlfXn", 
%4 .101f\n", 
%4 . lOlf\n", 
%4.lOlf\n", 
%4 . lOlf\n", 
%4 . lOlf\n", 
%4 . lOlf\n", 
%4 . lOlf\n", 
%4 . lOlf\n", 
%4.lOlfXn", 
%4 . lOlf\n", 
%4 . lOlf\n", 
%4 .lOlfXn", 
%4 . lOlf\n", 
%4 .lOlfXn", 
%4 . lOlf\n“, 
%4.lOlfXn", 
%4 .101f\n", 
%4 . lOlf\n", 
%4 .101f\n", 
%4 .101f\n", 
%4.101f\n", 
%4.101f\n", 
%4.101f\n", 
%4.101f\n", 
%4 . 101f\n", 
%4 . lOlfXn", 


b_in); 

b_wg_constant); 
b_wg1); 
b_wg2); 

b_tooth_tip_constant); 

b_tooth_tipl); 

b_tooth_tip2)/ 

cyclic_friction_amplitude); 

cyclic_friction_phase); 

stiction_torque); 

stiction_vel); 

b_out); 

mu_save); 

tooth_angle); 

N) ; 

cyclic_amplitude); 
cyclic_phase); 
error_amplitudeO); 
error_amplitudel) ; 
error_amplitude2) ; 
error_phaseO); 
error_phasel); 
error_phase2); 
motor_kt); 
motor_kb); 
amp_resistance) ; 
irequested); 


%4.lOlfXn", xstart_sav[POS_IN_INDEX ) ) ; 
%4.101f\n", xstart_sav[POS_OUT_INDEX ] ) ; 
%4.101f\n", xstart_sav[VEL_IN_INDEX]); 
%4.101f\n", xstart_sav[VEL_OUT_INDEX]); 

%4.lOlfXn", initial_time); 

%4.101f\n", final_time); 

%4.lOlfXn", step)/ 


/* Function print_variables () prints the names and values of some 
* important torques, velocities, and positions. */ 
void print variables() 


printf("\n"); 
printf("\n pos_in: 
printf("\n pos_out: 
printf("\n pos_wg: 
printf("\n pos_fs: 
printf("\n pos_cs: 
printf("\n pos_n_wg 
printf("\n"); 
printf("\n vel_in: 
printf("\n vel_out: 
printf("\n vel_wg: 
printf ("\n vel_fs: 
printf("\n vel_cs: 
printf ("\n"); 
printf("\n torque_motor 
printf("\n torque_b_in: 
printf("Xn torque_b_out: 
printf("\n torque_wg: 
printf("\n torque fs: 


%20.151f", v(POS_IN)); 
%20.151f", v[POS_OUT]); 

%20.151f", v[POS_WG]) ; 
%20.151f", v[POS_FS]); 
%20.151f", v1POS_CS]); 
%20.151f", v[POS_N_WG]); 

%20.151f", v [ V E L__ IN 1 ) ; 
%20.151f", v(VEL_OUT]); 
%20.151f", v[VEL_WG]); 
%20.151f", v[VEL_FS]); 
%20.151f", v[VEL_CS]); 

%20.151f", v[TORQUE_MOTOR]); 
%20.151f", v[TORQUE_B_IN]) ; 
%20.151f", v[TORQUE_B_OUT]); 
%20.151f", v[TORQUE_WG]); 
%20.151f", v[TORQUE FS]); 
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printf("\n torque_cs: 
printf ("\n torque_k: 
printf("\n torque_cyclic: 
printf("\n tooth_surface_friction: 
printf ("\n tooth_surf_total_loss: 
printf ("\n tooth_surface_normal: 
printf("\n torque_hd_friction: 
printf ("\n"); 
printf("\n current time: 
printf ("\n"); 
fflush(stdout); 


%20.151f ", v[TORQUE_CS]) ; 

%20.151f", v[TORQUE_K]); 

%20.151f", v[TORQUE_CYCLIC]) ; 

%20.151f", v[TORQUE_TOOTH_SURFACE_FRICTION]); 

%2 0.151f ", v[TORQUE_TOOTH_SURFACE_TOTAL_LOSS]) ; 
%20.151f", v[TORQUE_TOOTH_SURFACE NORMAL]); 

%20.151f", v[TORQUE_HD_FRICTION]); 

%20.151f", current time); 


/* Function export_vectors() takes a file pointer, two pointers 

* to vectors of 'kount' doubles and writes them to a 

file, each line containing an x and y coordinate. If the slow- 
plotting option is selected, the filename will not be printed 

* to the first line of the file. */ 

void export_vectors(outfile, xvector, yvector, points, filename) 

FILE *outfile; 

double *xvector; 

double *yvector; 

int points; 

char filename[]; 

{ 

int i; 

printf ("\n Exporting vector to file...."); 
if(!fast_plot) 

fprintf(outfile, "\" %s\n", filename); 
for(i = 1; i <= points; i++) 

fprintf(outfile, "%12.101f %12.lOlf\n", xvector[i], yvector[i]); 
printf ("done.\n"); 


/* Function generate_gnuplot_file() is called when the fast_plot option 
is selected and creates a command file in the data directory which 
is executed in gnuplot by a system command. */ 
void generate_gnuplot_file(filename) 
char filename]]; 

{ 

int i; 

FILE *gnufile; 

create_file(4gnu file, outdatapath, filename); 

fprintf(gnufile, "cd '%s'\n", outdatapath); 
fprintf (gnufile, "set xlabel 'time (sec) '\n") ; 

for(i =0; i < num_selected output var; i++) 

{ 

fprintf(gnufile, "set title '%s'\n", graph_title[i]); 
fprintf(gnufile, "set ylabel '%s'\n", graph_units[i]); 
fprintf(gnu file, "plot '%s'\n", output_datafile[i]); 
fprintf(gnufile, "pause -l\n"); 

} 

close_file(gnufile) ; 
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/* This is file hd model.h. 


* It is the header file for program hd model.c. */ 

/* Define a bunch of constants. */ 

♦define SHOULDER 

0 

♦define ELBOW 

1 

♦define WRIST 

2 

♦define JOINT TO TEST 

WRIST 

♦define NUM OF AXIS 

3 

♦define STIFFNESS POINTS 

1000 

♦define TRUE 

(1) 

♦define FALSE 

(0) 

♦define TINY_VEL 

(1.0e-20) 

/* Constants specifying the 

position of the local variables 

* in the state vector. */ 


♦define POS IN INDEX 

1 

♦define POS OUT INDEX 

2 

♦define VEL IN INDEX 

3 

♦define VEL OUT INDEX 

4 

♦define ENERGY 

5 

/* The order of the system 

(number of state variables). */ 

♦define ORDER 

4 

/* Define index numbers for 

two different harmonic-drive models. */ 

♦define ROTARY 

0 

♦define GEAR TOOTH 

1 

♦define HD_MODELS 

2 

/* Define the degrees-to radians conversion factor used in the equations. */ 

♦define PI 

3.14159265359 

♦define CONV 

(PI/180.0) 

/ Define mathematical functions to translate Mathematica's C-formatted 

* equations into a form understandable to the compiler. */ 

♦define Power(x,y) 

pow((double) x, (double) y) 

♦define Sin(x) 

sin((double) x) 

♦define Cos(x) 

cos((double) x) 

/* Define the simulation accuracy variable. */ 

♦define EPS 

0.00001 

/* Define a arbitrary intial 

step size and a minimum step size. */ 

♦define INITIAL STEP 

0.0001 

♦define MINIMUM STEP 

0.0000000000000000001 

/* Define constants describi 

ng the amplifier voltage-current curve */ 

♦define IMAX 

20.0 

♦define VMAX 

o 

o 

♦define VI SLOPE 

-1.0 

♦define AMP DAMPING 

0.5 

/* Define the number of data 

runs to be made if the -plot option is selected */ 

1ffdetine SHLD_NUM OF DATA RUNS 10 

|♦define ELB NUM OF DATA RUNS 

8 

♦define WRIST_NUM_OF_DATA RUNS 8 

/* Define the indicies for the vector to hold all variables that can be 

* selected for output data plots. */ 
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am 


#define 

POS WG 

0 

#define 

POS FS 

1 

#define 

POS CS 

2 

#define 

POS TOOTH BASE 

3 

#define 

POS ERR 

4 

#define 

POS TOOTH TIP 

5 

#define 

POS WG SURFACE 

6 

♦define 

POS TOOTH SURFACE 

7 

♦define 

POS WG RELATIVE 

8 

♦define 

POS K 

9 

♦define 

POS N WG 

10 

♦define 

POS N FS 

11 

♦define 

POS N CS 

12 

♦define 

VEL WG 

13 

♦define 

VEL FS 

14 

♦define 

VEL CS 

15 

♦define 

VEL TOOTH BASE 

16 

♦define 

VEL ERR 

17 

♦define 

VEL TOOTH TIP 

18 

♦define 

VEL_WG SURFACE 

19 

♦define 

VEL TOOTH SURFACE 

20 

♦define 

VEL WG RELATIVE 

21 

♦define 

VEL K 

22 

♦define 

VEL N WG 

23 

♦define 

VEL N FS 

24 

♦define 

VEL N CS 

25 

♦define 

VEL HD FRICTION 

26 

♦define 

TORQUE WG 

27 

♦define 

TORQUE FS 

28 

♦define 

TORQUE_CS 

29 

♦define 

TORQUE_INPUT NORMAL 

30 

♦define 

TORQUE_CYCLIC 

31 

♦define 

TORQUE_WG NORMAL 

32 

♦define 

TORQUE WG FRICTION 

33 

♦define 

TORQUE ERR IN 

34 

♦define 

TORQUE_K 

35 

♦define 

TORQUE_TOOTH TIP NORMAL 

36 

♦define 

TORQUE_TOOTH TIP FRICTION 

37 

♦define 

TORQUE_TOOTH_SURFACE NORMAL 

38 

♦define 

TORQUE_TOOTH_SURFACE TOTAL LOSS 

39 

♦define 

TORQUE_TOOTH SURFACE COULOMB 

40 

♦define 

TORQUE_TOOTH SURFACE FRICTION 

41 

♦define 

TORQUE CYCLIC FRICTION 

42 

♦define 

TORQUE OUTPUT NORMAL 

43 

♦define 

TORQUE N WG 

44 

♦define 

TORQUE N FS 

45 

♦define 

TORQUE_N CS 

46 

♦define 

TORQUE_HD_FRICTION 

47 

/* The positions, velocities, and torques 

outside the harmonic drive. */ 

♦define 

POS IN 

48 

♦define 

POS OUT 

4 9 

♦define 

VEL IN 

50 

♦define 

VEL_OUT 

51 

♦define 

TORQUE MOTOR 

52 

♦define 

TORQUE_B IN 

53 

♦define 

TORQUE B OUT 

54 

/* Sensor variables. */ 

♦define INPUT_POSITION SENSOR 

55 

♦define 

OUTPUT POSITION SENSOR 

56 

♦define 

INPUT_VELOCITY SENSOR 

57 

♦define 

OUTPUT_VELOCITY SENSOR 

58 

♦define 

CURRENT SENSOR 

59 
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#define TORQUE SENSOR 

60 


♦define POSITION ERROR 

61 


/* Parameter functions which can be plotted. */ 
♦define K1 

62 


♦define K2 

63 


♦define ERFN 

64 


♦define DERF 

65 


/* Energy Variables. */ 

♦define TOTAL ENERGY 

66 


♦define TOTAL_KINETIC ENERGY 

67 


♦define TOTAL_POTENTIAL ENERGY 

68 


♦define TOTAL INPUT ENERGY 

69 


♦define TOTAL LOST ENERGY 

70 


♦define MOTOR ENERGY 

71 


♦define KINETIC ENERGY IN 

72 


♦define KINETIC ENERGY OUT 

73 


♦define SPRING ENERGY 

74 


♦define CYCLIC ENERGY 

75 


♦define INPUT DAMPING ENERGY 

76 


♦define WG_FRICTION ENERGY 

77 


♦define TOOTH^TIP FRICTION ENERGY 

78 


♦define TOOTH_SURFACE TOTAL LOSS ENERGY 

79 


♦define TOOTH_SURFACE COULOMB ENERGY 

80 


♦define TOOTH_SURFACE FRICTION ENERGY 

81 


♦define HD_FRICTION ENERGY 

82 


♦define OUTPUT_DAMPING ENERGY 

83 


/* Define the number of variables in the output data vector. */ 


♦define NUM_OF OUTPUT VAR 

84 


/* Declare the variables for the input parameters 

used in the harmonic drive 

* equations. */ 

double kl_output_constant, kl_output amplitude, kl 

output phase; 


double k2_output_constant, k2 output amplitude, k2 

output phase; 


double hysteresis width; 



double error_amplitudeO, error phaseO; 
double error_amplitudel, error phasel; 



double error_amplitude2, error phase2; 
double cyclic_amplitude, cyclic phase; 



double b_wg_constant, b wgl, b wg2; 
double b_tooth_tip_constant, b_tooth_tipl, b tooth 
double cyclic_friction amplitude, cyclic friction 
double mu, mu save; 

_tip2; 

phase; 


double b_totall, b total2; 
double tooth_angle, N; 



double stiction_torque, stiction vel; 



/* Declare the variables for the joint parameters, 
double inertia_in, inertia out; 
double b_in, b out; 

*/ 


double motor_kt, motor kb, amp resistance; 



/* Requested step in motor current. */ 
double irequested; 



/* and final time, the stepsize interval to save results 

* / 

double initial_time, final time, step; 


/* Define the initial-condition vectors.*/ 
double xstart[ORDER+1] , xstart sav[ORDER + 1]; 




for 
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* each simulation time step. */ 
double v[NUM_OF_OUTPUT_VAR]; 

/* Define output vectors to be allocated dynamically. */ 
double **output_data, *time; 

/* Define the parameters that are calculated from the input parameters 

* for the rotary harmonic drive model. */ 
double kl_constant, kl_amplitude, kl_phase; 
double k2_constant, k2_amplitude, k2_phase; 
double stiction_torque_hd, stiction_vel_hd; 

double cyclic_friction_amp_hd, cyclic_friction_phase hd; 
double b_hd_constant, b_hdl, b_hd2; 

/* Define the parameters that are calculated from the input parameters 

* for the tooth-rubbing harmonic drive model. */ 
double tanl, cosl, sinl, tan2, cos2, sin2; 
double wg angle; 

double stiction_torque_tooth, stiction_vel_tooth; 

double cyclic_friction_amp_tooth, cyclic_friction_phase_tooth; 

double b_tooth_surface_constant, b_tooth_surfacel, b_tooth_surface2; 
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For most of the operating range of the harmonic-drive testing equipment, the 
amplifiers which supply current to the DC motors can sustain current values up to the their 
maximum rated current. However, at high motor velocities, back EMF generated by the 
motor can saturate the amplifier voltage which results in a reduction of the maximum 
current the amplifier can deliver. This saturation behavior can be captured by developing a 
few simple equations which incorporate the velocity and current limits of the amplifiers 
with some simplified dynamics of the motor-amplifier circuit. 


Identical Aerotech model 4020-LS linear current-amplifiers were used for the three 
harmonic-drive testing stations. The voltage and current operating range of these amplifiers 
is shown in figure C.l. Using this plot, the upper and lower voltage limits of the amplifier 
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can be calculated for any desired current within the given current range, as illustrated in 
figure C.2. Given this upper and lower voltage-bound, an upper and lower current-bound 
can also be found by considering the simplified amplifier circuit shown in figure C.3. 
Summing the voltage around the circuit yields the equation 


Vamp ~ V motor + Vr , Or 


(C.l) 


Vamp ~ KbO) motor 


lamp R 


(C.2) 


where K b is the back EMF constant of the motor, co motor is the motor velocity, and R is the 
resistance across the motor terminals. Resistance values which provided good results are 
listed in Table C.l for the three harmonic-dnve testing stations. By solving equation C.2 
for the amplifier current, i a mp* and substituting the voltage upper and lower-bounds as 
determined in figure C.2, an upper and lower bound for the amplifier current can be 
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derived: 


1 


amp max 


Vamp max " Kb^motor 


R 


, and 


; _ ^amp min ' ^b^motor 

*amp min-- 


(C.3) 

(C.4) 


Therefore, for a given motor operating speed and a desired current value, the instantaneous 
current limit on the amplifiers can be determined simply by (1) using the amplifier operating 
specifications to calculate the voltage limits for the requested current, and (2) substituting 
these voltage limits into equations (C.3) and (C.4). This model was incorporated 
successfully into a numerical simulation which calculated the amplifier current on every 
time-step. 
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Appendix D: Experimental Time- 

Response Data Plots 


This appendix contains a complete set of all of the dynamic-response plots collected 
from the three harmonic-drive testing stations. For each harmonic-drive, input and output 
position were measured and used to derive input and output velocity data. The output 
torque on each transmission was monitored by the torque sensors while the input torque 
was determined from the motor-current sensor. Additionally, from the input and output 
position measurements, the dynamic position-error across the transmission was calculated. 
A summary of all of the plots shown in this appendix and their corresponding plot numbers 
is given in table D.l. For each of the plots listed in this table, except the dynamic position- 
error figure, time-response curves are plotted for several different motor-current step- 
commands. Depending on the visual appearance of the data, the different response curves 
on each of these plots are either plotted on a single axis or plotted on separate axes. For the 
sake of clarity, the output-velocity response data for each harmonic drive is plotted in both 
of these ways. In the case of the dynamic position-error plots, a single time-response 
curve is shown, rather than several, in order to generalize about behavior for different step 
inputs. A complete discussion of all of these plots can be found in section 2.6 of this 
document. 
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Table D.1: Summary of experimental time-response data-plot numbers 


. 


input Velocity 

.snsns^sv.v.v^vavaw.w.v.v.v.v.-:?:-: 


HHMi 


OuiptigTo rque |j|| 


Dynamic Position Error 


Joint 

llll 

Joint 2 

Joint 3 

D.1.1 

D.2.1 

D.3.1 

D.1.2 

D.2.2 

D.3.2 

D.1.3 

D.2.3 

D.3.3 

D.1.4 

D.2.4 

D.3.4 

D.1.5 

D.2.5 

D.3.5 

D.1.6 

D.2.6 

D.3.6 

D.1.7 

D.2.7 

D.3.7 

D.1.8 
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Appendix D: Experimental Data Plots for Joint 1 
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Appendix D: Experimental Data Plots for Joint 2 
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input velocity (degrees/second) 
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Input Motor Current (amperes) 
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output velocity (degrees/second) 
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Figure D.2.8: Joint 2 dynamic position error at 2.6 amps resonance 
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Appendix D: Experimental Data Plots for Joint 3 
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input velocity (degrees/second) 
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Figure D.3.4: Joint 3 output-position step-responses 
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output velocity (degrees/second) 
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Figure D.3.5: Joint 3 output-velocity step-responses 
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time (seconds) 


Figure D.3.8: Joint 3 dynamic position error at 0.48 amps resonance 
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Appendix E: Simulated Time- 

Response Plots for the 
Harmonic-Drive Model 
with Friction, Compliance 
and Kinematic Error 


This appendix contains a complete collection of all of the dynamic-response plots 
generated for each testing station using the harmonic-drive model with friction, compliance, 
and kinematic error. This model is derived and these results are discussed in section 3.6 of 
this document. The format of these illustrations is identical to the experimental time- 
response results presented in appendix D so that direct comparisons can be made easily. 
Using the numerical simulation described in section 3.1, values for the position, velocity, 
and torque in the model were calculated, and from the input and output position data, the 
dynamic position-error across the transmission was derived. A summary of all of the plots 
shown in this appendix and their corresponding plot numbers is given in table E. 1. For 
each of the plots listed in this table except the dynamic position-error figure, time-response 
curves are plotted for several different motor-current step-commands. Depending on the 
visual appearance of the data, the different response curves on each of these plots are either 
plotted on a single axis or plotted on separate axes. For the sake of clarity, the output- 
velocity response data for each harmonic drive is plotted in both of these ways. In the case 
of the dynamic position-error plots, a single time-response curve is shown, rather than 
several, from which generalizations about the behavior of dynamic position-error in all 
operating ranges are made. 
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Table E.1: Summary of simulated time-response plot numbers 
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output position (degrees) 



Figure E.1.4: Joint 1 output-position simulated step-responses 
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Figure E.1.5: Joint 1 output-velocity simulated step-responses 
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Appendix E: Model 4 Simulated Plots for Joint 1 
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Figure E.1.7: Joint 1 output-torque simulated step-responses 
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Figure E.1.8: Joint 1 simulated dynamic position error at 4.0 amps 
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Figure E.2.4: Joint 2 output-position simulated step-responses 
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Figure E.2.5: Joint 2 output-velocity simulated step-responses 
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Figure E.2.7: Joint 2 output-torque simulated step-responses 
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Figure E.2.8: 


simulated dynamic position error at 2.2 amps 
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output position (degrees) 



Figure E.3.4: Joint 3 output-position simulated step-responses 
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Figure E.3.8: Joint 3 simulated dynamic position error at 0.44 amps 
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Appendix F: Simulated Time- 

Response Plots for the 
Harmonic-Drive Model 
Including Gear-Tooth 
Geometry 


This appendix contains a complete collection of all of the dynamic-response plots 
generated for each testing station using the harmonic-drive model including gear-tooth 
geometry. This model is derived and these results are discussed in section 3.7 of this 
document. The format of these illustrations is identical to the experimental time-response 
results presented in appendix D so that direct comparisons can be made easily. Using the 
numerical simulation described in section 3.1, values for the position, velocity, and torque 
in the model were calculated, and from the input and output position data, the dynamic 
position-error across the transmission was derived. A summary of all of the plots shown 
in this appendix and their corresponding plot numbers is given in table F. 1. For each of 
the plots listed in this table except the dynamic position-error figure, time-response curves 
are plotted for several different motor-current step-commands. Depending on the visual 
appearance of the data, the different response curves on each of these plots are either 
plotted on a single axis or plotted on separate axes. For the sake of clarity, the output- 
velocity response data for each harmonic drive is plotted in both of these ways. In the case 
of the dynamic position-error plots, a single time-response curve is shown, rather than 
several, from which generalizations about the behavior of dynamic position-error in all 
operating ranges are made. 
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Figure F.1.4: Joint 1 output position simulated step-responses 
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Figure F.1.5: Joint 1 output-velocity simulated step-responses 
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Figure F.1.7: Joint 1 output-torque simulated step-responses 
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Figure F.2.3: Joint 2 input-current simulated step-responses 
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Figure F.2.4: Joint 2 output-position simulated step-responses 
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Figure F.2.7: Joint 2 output-torque simulated step-responses 
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Figure F.2.8: Joint 2 simulated dynamic position-error at 2.6 amps 
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Figure F.3.3: Joint 3 input-current simulated step-responses 
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ix F: Model 5 Simulated Plots for Joint 3 
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Appendix F: Model 5 Simulated Plots for Joint 3 
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endix F: Model 5 Simulated Plots for Joint 3 



Figure F.3.8: Joint 3 simulated dynamic position error at 0.44 amps 
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